solvingoperations.c

cfd求解器使用与gmsh网格的求解

#define RCSID "$Id: SolvingOperations.c,v 1.79 2006/03/13 22:48:21 geuzaine Exp $"
/*
 * Copyright (C) 1997-2006 P. Dular, C. Geuzaine
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA.
 *
 * Please report all bugs and problems to <[e-mail]>.
 *
 * Contributor(s):
 *   David Colignon
 *   Johan Gyselinck
 *   Ruth Sabariego
 */

#include "GetDP.h"
#include "Treatment_Formulation.h"
#include "GeoData.h"
#include "DofData.h"
#include "Parser.h"
#include "Init_Problem.h"
#include "Cal_Quantity.h"
#include "Tools.h"
#include "Data_DefineE.h"
#include "Numeric.h"
#include "Get_DofOfElement.h"
#include "CurrentData.h"
#include "ExtendedGroup.h"
#include "MovingBand2D.h"
#include "Magic.h"
#include "F_FMM.h"
#include "F_FMMOperations.h"
#include "F_FMM_DTA.h"
#include "BF_Function.h"

static int  Flag_IterativeLoop = 0 ;  /* Attention: phase de test */
static int  Flag_NextThetaFixed = 0 ;  /* Attention: phase de test */
static int  Init_Update = 0 ; /* provisoire */

struct Group * Generate_Group = NULL;

/* Johan: it would be nice to get rid opf these globals */
int Flag_RHS = 0, *DummyDof ;
double **MH_Moving_Matrix = NULL ; 
int MH_Moving_Matrix_simple = 0 ;
int MH_Moving_Matrix_probe = 0 ;
int MH_Moving_Matrix_separate = 0;
Tree_T  * DofTree_MH_moving ;

/*
static FILE * FilePWM ;
*/

/* ------------------------------------------------------------------------ */
/*  I n i t _ O p e r a t i o n O n S y s t e m                             */
/* ------------------------------------------------------------------------ */

void  Init_OperationOnSystem(char                * Name,
			     struct Resolution   * Resolution_P,
			     struct Operation    * Operation_P ,
			     struct DofData      * DofData_P0,
			     struct GeoData      * GeoData_P0,
			     struct DefineSystem ** DefineSystem_P,
			     struct DofData      ** DofData_P,  
			     struct Resolution   * Resolution2_P) {

  int i ;

  GetDP_Begin("Init_OperationOnSystem");

  *DefineSystem_P = (struct DefineSystem*)
    List_Pointer(Resolution_P->DefineSystem,Operation_P->DefineSystemIndex) ;
  *DofData_P = DofData_P0 + Operation_P->DefineSystemIndex ;
  Dof_SetCurrentDofData(Current.DofData = *DofData_P) ;
  Current.NbrHar = Current.DofData->NbrHar ;

  Geo_SetCurrentGeoData(Current.GeoData = GeoData_P0 + (*DofData_P)->GeoDataIndex) ;

  if((*DefineSystem_P)->DestinationSystemName && 
     (*DefineSystem_P)->DestinationSystemIndex == -1){

    if(Resolution2_P){ /* pre-resolution */
      if ((i = List_ISearchSeq(Resolution2_P->DefineSystem, 
			       (*DefineSystem_P)->DestinationSystemName,
			       fcmp_DefineSystem_Name)) < 0)
	Msg(GERROR, "Unknown DestinationSystem (%s) in System (%s)", 
	    (*DefineSystem_P)->DestinationSystemName, (*DefineSystem_P)->Name) ;
      (*DefineSystem_P)->DestinationSystemIndex = i ;      
      Dof_DefineUnknownDofFromSolveOrInitDof(DofData_P) ;
    }
    else { /* a changer !!! */
      if ((i = List_ISearchSeq(Resolution_P->DefineSystem, 
			       (*DefineSystem_P)->DestinationSystemName,
			       fcmp_DefineSystem_Name)) < 0)
	Msg(GERROR, "Unknown DestinationSystem (%s) in System (%s)", 
	    (*DefineSystem_P)->DestinationSystemName, (*DefineSystem_P)->Name) ;
      (*DefineSystem_P)->DestinationSystemIndex = i ;      
    }
  }
  
  Msg(OPERATION, "%s[%s]",
      Name?Name:Get_StringForDefine(Operation_Type, Operation_P->Type),
      (*DefineSystem_P)->Name) ;

  GetDP_End ;
}


/* ------------------------------------------------------------------------ */
/*  I n i t _ S y s t e m D a t a                                           */
/* ------------------------------------------------------------------------ */

void  Init_SystemData(struct DofData * DofData_P, int Flag_Jac) {

  GetDP_Begin("Init_SystemData");

  if (DofData_P->Flag_Init[0] < 1) {
    DofData_P->Flag_Init[0] = 1 ;
    LinAlg_CreateSolver(&DofData_P->Solver, DofData_P->SolverDataFileName) ;
    LinAlg_CreateMatrix(&DofData_P->A, &DofData_P->Solver,
			DofData_P->NbrDof, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part, DofData_P->Nnz) ;
    LinAlg_CreateVector(&DofData_P->b, &DofData_P->Solver, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part) ;
  }

  /* GenerateOnly: Taking advantage of the invariant parts of the matrix in every time-step */

  if(DofData_P->Flag_InitOnly[0]==1){
    DofData_P->Flag_InitOnly[0] = 2;    
    Msg(INFO,"Initializing System {A1,b1}");
    LinAlg_CreateMatrix(&DofData_P->A1, &DofData_P->Solver,
			DofData_P->NbrDof, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part, DofData_P->Nnz) ;
    LinAlg_CreateVector(&DofData_P->b1, &DofData_P->Solver, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part) ;
  }
  
  if(DofData_P->Flag_InitOnly[1]==1){
    DofData_P->Flag_InitOnly[1] = 2;    
    Msg(INFO,"Initializing System {A2,b2}");
    LinAlg_CreateMatrix(&DofData_P->A2, &DofData_P->Solver,
			DofData_P->NbrDof, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part, DofData_P->Nnz) ;
    LinAlg_CreateVector(&DofData_P->b2, &DofData_P->Solver, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part) ;
  }
  
  if(DofData_P->Flag_InitOnly[2]==1){
    DofData_P->Flag_InitOnly[2] = 2;       
    Msg(INFO,"Initializing System {A2,b2}");
    LinAlg_CreateMatrix(&DofData_P->A3, &DofData_P->Solver,
			DofData_P->NbrDof, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part, DofData_P->Nnz) ;
    LinAlg_CreateVector(&DofData_P->b3, &DofData_P->Solver, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part) ;
  }

  if (DofData_P->Flag_Init[0] < 2 && Flag_Jac) {
    DofData_P->Flag_Init[0] = 2 ;
    LinAlg_CreateMatrix(&DofData_P->Jac, &DofData_P->Solver,
			DofData_P->NbrDof, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part, DofData_P->Nnz) ;
    LinAlg_CreateVector(&DofData_P->res, &DofData_P->Solver, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part) ;
    LinAlg_CreateVector(&DofData_P->dx, &DofData_P->Solver, DofData_P->NbrDof,
			DofData_P->NbrPart, DofData_P->Part) ;
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  C a l _ S o l u t i o n E r r o r                                       */
/* ------------------------------------------------------------------------ */

void  Cal_SolutionError(gVector *dx, gVector *x, int diff, double *MeanError) {
  int     i, n;
  double  valx, valdx, errsqr = 0., xmoy = 0., dxmoy = 0., tol ;

  GetDP_Begin("Cal_SolutionError");

  if(gSCALAR_SIZE == 2)
    Msg(WARNING, "FIXME: Cal_SolutionError might return strange results"
	" in complex arithmetic");
  
  LinAlg_GetVectorSize(dx, &n);

  for (i=0 ; i<n ; i++) {
    LinAlg_GetAbsDoubleInVector(&valx, x, i) ; 
    LinAlg_GetAbsDoubleInVector(&valdx, dx, i) ; 
    xmoy += valx ;
    if(diff) dxmoy += (valdx-valx) ;
    else     dxmoy += valdx ;
  }
  xmoy  /= (double)n ;
  dxmoy /= (double)n ;

  if (xmoy > 1.e-30) {
    tol = xmoy*1.e-10 ;
    for (i=0 ; i<n ; i++){
      LinAlg_GetAbsDoubleInVector(&valx, x, i) ;
      LinAlg_GetAbsDoubleInVector(&valdx, dx, i) ;
      if(diff){
	if (valx > tol) errsqr += fabs(valdx-valx)/valx ;
	else 	        errsqr += fabs(valdx-valx) ;
      }
      else{
	if (valx > tol) errsqr += valdx/valx ;
	else 	        errsqr += valdx ;
      }
    }
    *MeanError = errsqr/(double)n ;
  }
  else{
    if (dxmoy > 1.e-30) 
      *MeanError = 1. ;
    else
      *MeanError = 0. ;
  }

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  C a l _ S o l u t i o n E r r o r X                                     */
/* ------------------------------------------------------------------------ */

void  Cal_SolutionErrorX(int Nbr, double * xNew, double * x, double * MeanError) {
  int i;
  double errsqr = 0., xmoy = 0., dxmoy = 0., tol ;

  GetDP_Begin("Cal_SolutionErrorX");

  if(gSCALAR_SIZE == 2)
    Msg(GERROR, "FIXME: Cal_SolutionErrorX might return strange results"
	" in complex arithmetic");

  for (i = 0 ; i < Nbr ; i++) {
    xmoy  += fabs( x[i])/(double)Nbr ;
    dxmoy += fabs(xNew[i]-x[i])/(double)Nbr ;
  }

  if (xmoy > 1.e-30) {
    tol = xmoy*1.e-10 ;

    for (i = 0 ; i < Nbr ; i++)
      if ( fabs(x[i]) > tol ) errsqr += fabs((xNew[i]-x[i]) / x[i]) ;
      else                    errsqr += fabs(xNew[i]-x[i]) ;

    *MeanError = errsqr / (double)Nbr ;
  }
  else
    if (dxmoy > 1.e-30)  *MeanError = 1. ;
    else                 *MeanError = 0. ;

  GetDP_End ;
}

/* ------------------------------------------------------------------------ */
/*  T r e a t m e n t _ O p e r a t i o n                                   */
/* ------------------------------------------------------------------------ */

void  Treatment_Operation(struct Resolution  * Resolution_P,
                          List_T             * Operation_L, 
                          struct DofData     * DofData_P0,
                          struct GeoData     * GeoData_P0,
                          struct Resolution  * Resolution2_P,
                          struct DofData     * DofData2_P0) {

  int     i, j, k, l ;
  double  d, d1, d2, *Scales ;
  int     Nbr_Operation, Nbr_Sol, i_Operation, Num_Iteration ;
  int     Flag_Jac, Flag_CPU, Flag_Binary = 0, Flag_FMMDA = 0 ;
  int     Save_TypeTime ;
  double  Save_Time, Save_TimeImag, Save_DTime, Save_TimeStep ;
  double  MeanError, RelFactor_Modified ;
  char    *str;
  char    ResName[MAX_FILE_NAME_LENGTH], ResNum[MAX_STRING_LENGTH] ;
  char    FileName[MAX_FILE_NAME_LENGTH], FileFMM[MAX_FILE_NAME_LENGTH];
  char    FileName_exMH[MAX_FILE_NAME_LENGTH];
  gScalar tmp ;

  struct Operation      * Operation_P ;
  struct DefineSystem   * DefineSystem_P ;
  struct DofData        * DofData_P, * DofData2_P ;
  struct Solution       * Solution_P, Solution_S ;
  struct PostOperation  * PostOperation_P ;
  struct PostProcessing * PostProcessing_P ;
  struct Dof            Dof, * Dof_P ;
  struct Value          Value, far ;

  double  **DTA ;
  FILE * MomFMM ;
  int iDof, iEqu, N ;

  static int RES0 = -1 ;

  /* adaptive relaxation */
  gVector x_Save;
  int NbrSteps_relax;
  double  Norm;
  double Frelax, Frelax_Opt, Error_Prev;
  int istep;

  int Nbr_Formulation, Index_Formulation ;
  struct Formulation * Formulation_P ; 

  int iTime ;
  double *Val_Pulsation ; 
  double hop[NBR_MAX_HARMONIC][NBR_MAX_HARMONIC] ;
  double DCfactor ;

  int NbrHar1, NbrHar2, NbrDof1, NbrDof2 ;
  double dd ;
  int NumDof, iMat ;
  int NbrFMMEqu;
  struct FMMmat *FMMmat_P0 ;

  int row_old, row_new, col_old, col_new;

  double aii, ajj;
  int nnz__;

  List_T  * DofList_MH_moving ;
  static int NbrDof_MH_moving;
  static int *NumDof_MH_moving ; 
  static struct Dof ** Dof_MH_moving ; 

  GetDP_Begin("Treatment_Operation");

  Nbr_Operation = List_Nbr(Operation_L) ;

  for (i_Operation = 0 ; i_Operation < Nbr_Operation ; i_Operation++) {
    Operation_P = (struct Operation*)List_Pointer(Operation_L, i_Operation) ;

    Flag_CPU = 0 ;  
    Flag_Jac = 0 ;

    switch (Operation_P->Type) {

      /*  -->  S y s t e m C o m m a n d              */
      /*  ------------------------------------------  */