solvingoperations.c

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

    case OPERATION_SYSTEMCOMMAND :
      system(Operation_P->Case.SystemCommand.String);
      break ;

      /*  -->  G e n e r a t e                        */
      /*  ------------------------------------------  */

    case OPERATION_GENERATEJAC :  Flag_Jac  = 1 ;   
    case OPERATION_GENERATE :
      Init_OperationOnSystem(Get_StringForDefine(Operation_Type, Operation_P->Type),
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
			     &DefineSystem_P, &DofData_P, Resolution2_P) ;

      if ( Current.FMM.SystemIndex == Operation_P->DefineSystemIndex )
	if (Flag_FMM)
	  DefineSystem_P->Flag_FMM = Flag_FMM ;
	else 
	  Flag_FMM = DefineSystem_P->Flag_FMM ;
      else Flag_FMM = DefineSystem_P->Flag_FMM ;

      Current.TypeAssembly = ASSEMBLY_AGGREGATE ;
      
      Init_SystemData(DofData_P, Flag_Jac) ;

      if (Operation_P->Case.Generate.GroupIndex >= 0) 
	Generate_Group = (struct Group *) List_Pointer(Problem_S.Group, 
						       Operation_P->Case.Generate.GroupIndex) ;
      Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 0) ;

      if (Operation_P->Case.Generate.GroupIndex >= 0) Generate_Group = NULL ;

      Flag_CPU = 1 ;
      break ;

      /*  -->  G e n e r a t e S e p a r a t e        */
      /*  ------------------------------------------  */

    case OPERATION_GENERATESEPARATE :
      Init_OperationOnSystem("GenerateSeparate",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;

      if (Operation_P->Case.Generate.GroupIndex >= 0) 
	Generate_Group = (struct Group *) List_Pointer(Problem_S.Group, 
						       Operation_P->Case.Generate.GroupIndex) ;
      Current.TypeAssembly = ASSEMBLY_SEPARATE ;
      Init_Update = 0 ; /* modif... ! */

      Init_SystemData(DofData_P, Flag_Jac) ;
      Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 1) ;

      if (Operation_P->Case.Generate.GroupIndex >= 0) Generate_Group = NULL ;
      Flag_CPU = 1 ;
      break ;

      /*  -->  G e n e r a t e O n l y                */
      /*  ------------------------------------------  */

    case OPERATION_GENERATEONLYJAC :  Flag_Jac  = 1 ;
    case OPERATION_GENERATEONLY:
      Init_OperationOnSystem(Get_StringForDefine(Operation_Type, Operation_P->Type),
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
			     &DefineSystem_P, &DofData_P, Resolution2_P) ;

      if ( Current.FMM.SystemIndex == Operation_P->DefineSystemIndex )
	if (Flag_FMM)
	  DefineSystem_P->Flag_FMM = Flag_FMM ;
	else 
	  Flag_FMM = DefineSystem_P->Flag_FMM ;
      else Flag_FMM = DefineSystem_P->Flag_FMM ;

      if(DofData_P->Flag_Only < 2) DofData_P->Flag_Only += 1 ;
      
      DofData_P->OnlyTheseMatrices = Operation_P->Case.GenerateOnly.MatrixIndex_L ;
      
      if (DofData_P->Flag_Only <= 2)
	for (i = 0 ; i < List_Nbr(DofData_P->OnlyTheseMatrices); i++){
	  List_Read( DofData_P->OnlyTheseMatrices, i, &iMat);
	  switch(iMat){
	  case 1: DofData_P->Flag_InitOnly[0] = 1 ; break ;
	  case 2: DofData_P->Flag_InitOnly[1] = 1 ; break ;
	  case 3: DofData_P->Flag_InitOnly[2] = 1 ; break ;
	  }
	}

      Current.TypeAssembly = ASSEMBLY_AGGREGATE ;
      
      Init_SystemData(DofData_P, Flag_Jac) ;
      Generate_System(DefineSystem_P, DofData_P, DofData_P0, Flag_Jac, 0) ;
      Flag_CPU = 1 ;
      break;

      /*  -->  G e n e r a t e F M M G r o u p s      */
      /*  ------------------------------------------  */

    case OPERATION_GENERATEFMMGROUPS : Flag_FMM  = 1 ;       
      Init_OperationOnSystem(Get_StringForDefine(Operation_Type, Operation_P->Type),
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;
      
      DefineSystem_P->Flag_FMM = Flag_FMM ;      
      Current.FMM.SystemIndex = Operation_P->DefineSystemIndex ;
            
      if (List_Nbr(DefineSystem_P->FormulationIndex) == 1)
	List_Read(DefineSystem_P->FormulationIndex, 0, &Index_Formulation) ;
      else
	Msg(GERROR,"FMM not ready for working with more than 1 formulation per system!");

      Current.FMM.DivXYZIndex = Operation_P->Case.GenerateFMMGroups.DivXYZIndex ;
      Get_ValueOfExpressionByIndex(Operation_P->Case.GenerateFMMGroups.Dfar,
      				   NULL, 0., 0., 0., &far) ;
      /* far = Factor for determining the far groups */
      Current.FMM.far = far.Val[0] ;

      if (Operation_P->Case.GenerateFMMGroups.Precision == -1)
	Current.FMM.Precision = 0. ;
      else{
	Get_ValueOfExpressionByIndex(Operation_P->Case.GenerateFMMGroups.Precision,
				     NULL, 0., 0., 0., &Value) ; /* Precision */
	Current.FMM.Precision = Value.Val[0] ;
      }
     
      if (Operation_P->Case.GenerateFMMGroups.FlagDTA == -1)
	Flag_DTA = 0 ;
      else{
	Get_ValueOfExpressionByIndex(Operation_P->Case.GenerateFMMGroups.FlagDTA,
				     NULL, 0., 0., 0., &Value) ; /* DTA matrix, testing option */
	Flag_DTA = (int)Value.Val[0] ;
      }

      Current.FMM.NbrDir = 0 ;
      Get_InFMMGroupList( Index_Formulation, GeoData_P0+DofData_P->GeoDataIndex );
      
      strcpy(FileFMM, Name_Generic) ; strcat(FileFMM, "GrCen.pos") ;
      Geo_WriteFileFMMGroupsCenter(FileFMM) ;

      strcpy(FileFMM, Name_Generic);strcat(FileFMM, "Gr.pos") ;
      Geo_WriteFileFMMGroups( GeoData_P0+DofData_P->GeoDataIndex , FileFMM ) ;
      break;

      /*  -->  U p d a t e                            */
      /*  ------------------------------------------  */

    case OPERATION_UPDATE :
      Init_OperationOnSystem("Update",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;
      Update_System(DefineSystem_P, DofData_P, DofData_P0, 
		    Operation_P->Case.Update.ExpressionIndex) ;
      Flag_CPU = 1 ;
      break ;

      /*  -->  U p d a t e C o n s t r a i n t        */
      /*  ------------------------------------------  */

    case OPERATION_UPDATECONSTRAINT :
      Init_OperationOnSystem("UpdateConstraint",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;
      UpdateConstraint_System(DefineSystem_P, DofData_P, DofData_P0, 
			      Operation_P->Case.UpdateConstraint.GroupIndex,
			      Operation_P->Case.UpdateConstraint.Type) ;
      Flag_CPU = 1 ;
      break ;

      /*  -->  U p d a t e F M M                      */
      /*  ------------------------------------------  */

    case OPERATION_UPDATEFMMDATA : Flag_FMMDA = 1 ;
    case OPERATION_UPDATETRANSLATION :
      Init_OperationOnSystem(Get_StringForDefine(Operation_Type, Operation_P->Type),
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;

      ReSet_FMMGroupCenters();      

      ReGenerate_FMMGroupNeighbours(Flag_FMMDA);

      NbrFMMEqu = List_Nbr(Current.DofData->FMM_Matrix) ;
      FMMmat_P0 = (struct FMMmat*)List_Pointer(Current.DofData->FMM_Matrix, 0 ) ;

      for(iEqu = 0 ; iEqu < NbrFMMEqu ; iEqu ++ ){
	if((FMMmat_P0 + iEqu)->T != NULL){
	  Free((FMMmat_P0 + iEqu)->T);

	  (FMMmat_P0 + iEqu)->T = NULL;
	}
      } 

      GF_FMMTranslationValueAdaptive();
      /* GF_FMMTranslationValue(); */
      break;

      /*  -->  S o l v e                              */
      /*  ------------------------------------------  */

    case OPERATION_SOLVE :
      /*  Solve : A x = b  */
      Init_OperationOnSystem("Solve",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;

      if ( Current.FMM.SystemIndex == Operation_P->DefineSystemIndex )
	if (Flag_FMM) DefineSystem_P->Flag_FMM = Flag_FMM ;


      if (Flag_FMM && Flag_DTA){	
	LinAlg_GetMatrixSize(&DofData_P->A, &N, &N);
	FMM_DTAMatrix(N, &DTA);

	strcpy(FileFMM, Name_Generic) ;
	strcat(FileFMM, "DTA.mat") ;
	Print_FMM_MatrixDTA(N, DTA, FileFMM) ;
	
	for(iDof = 0 ; iDof < N ; iDof+=Current.NbrHar)
	  for(iEqu = 0; iEqu < N ; iEqu+=Current.NbrHar)
	    if (Current.NbrHar == 2)
	      LinAlg_AddComplexInMatrix(DTA[iDof][iEqu],DTA[iDof+1][iEqu+1],
					&DofData_P->A, iDof, iEqu, iDof+1, iEqu+1);
	    else   
	      LinAlg_AddDoubleInMatrix(DTA[iEqu][iDof], &DofData_P->A, iDof, iEqu );

	strcpy(FileFMM, Name_Generic) ;
	strcat(FileFMM, "FMM.mat") ;

	Msg(INFO,"Printing MoM + DTA matrix '%s'", FileFMM) ;
	MomFMM = fopen(FileFMM, "w" );
	LinAlg_PrintMatrix(MomFMM , &DofData_P->A) ;
	fclose(MomFMM);
      }
     
      if (DofData_P->Flag_Only){
	if(DofData_P->Flag_InitOnly[0]){
	  LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A1, &DofData_P->A);
	  LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b1, &DofData_P->b) ;
	}

	if(DofData_P->Flag_InitOnly[1]){
	  LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A2, &DofData_P->A) ;
	  LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b2, &DofData_P->b) ;
	}
	if(DofData_P->Flag_InitOnly[2]){
	  LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A3, &DofData_P->A) ;
	  LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b3, &DofData_P->b) ;	
	}

	LinAlg_AssembleMatrix(&DofData_P->A) ;
	LinAlg_AssembleVector(&DofData_P->b) ;
      }

      LinAlg_Solve(&DofData_P->A, &DofData_P->b, &DofData_P->Solver,
		   &DofData_P->CurrentSolution->x) ;

      Flag_CPU = 1 ;
      break ;

      /*  -->  S o l v e J a c                        */
      /*  ------------------------------------------  */

    case OPERATION_SOLVEJAC :
      /*  SolveJac : J(xn) dx = b(xn) - A(xn) xn ;  x = xn + dx  */
      Flag_Jac = 1 ;
      Init_OperationOnSystem("SolveJac",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;

      if ( Current.FMM.SystemIndex == Operation_P->DefineSystemIndex )
	if (Flag_FMM) DefineSystem_P->Flag_FMM = Flag_FMM ;

      if(DofData_P->Flag_Init[0] < 2)
	Msg(GERROR, "Jacobian system not initialized (missing GenerateJac?)");
      
      if (DofData_P->Flag_Only){
	if(DofData_P->Flag_InitOnly[0]){
	  LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A1, &DofData_P->A);
	  LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b1, &DofData_P->b) ;
	}

	if(DofData_P->Flag_InitOnly[1]){
	  LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A2, &DofData_P->A) ;
	  LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b2, &DofData_P->b) ;
	}
	if(DofData_P->Flag_InitOnly[2]){
	  LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A3, &DofData_P->A) ;
	  LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b3, &DofData_P->b) ;	
	}

	LinAlg_AssembleMatrix(&DofData_P->A) ;
	LinAlg_AssembleVector(&DofData_P->b) ;
      }

      LinAlg_AddMatrixMatrix(&DofData_P->Jac, &DofData_P->A, &DofData_P->Jac) ;
      LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, &DofData_P->res) ;
     
      if (Flag_FMM)
	LinAlg_FMMMatVectorProd(&DofData_P->CurrentSolution->x, &DofData_P->res) ;
      
      LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res) ;
      LinAlg_DummyVector(&DofData_P->res) ;
      LinAlg_Solve(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver, &DofData_P->dx) ;

      Cal_SolutionError(&DofData_P->dx, &DofData_P->CurrentSolution->x, 0, &MeanError) ;
      Msg(BIGINFO, "Mean error: %.3e  (after %d iteration%s)", 
	  MeanError, (int)Current.Iteration, ((int)Current.Iteration==1)?"":"s") ;

      Current.RelativeDifference += MeanError ;

      if (!Flag_IterativeLoop) {
        LinAlg_ProdVectorDouble(&DofData_P->dx, Current.RelaxationFactor, &DofData_P->dx) ;
      }
      else {  /* Attention: phase test ... Technique bricolee ... provisoire */
        if (Current.Iteration == 1. || MeanError < Current.RelativeDifferenceOld)
          LinAlg_ProdVectorDouble(&DofData_P->dx, Current.RelaxationFactor, &DofData_P->dx) ;
        else {
          RelFactor_Modified = Current.RelaxationFactor /
            (MeanError / Current.RelativeDifferenceOld) ;
          Msg(INFO, "RelFactor modified = %g", RelFactor_Modified) ;
          LinAlg_ProdVectorDouble(&DofData_P->dx, RelFactor_Modified, &DofData_P->dx) ;
          Cal_SolutionError(&DofData_P->dx, &DofData_P->CurrentSolution->x, 0, &MeanError) ;
          Msg(BIGINFO, "Mean error: %.3e", MeanError) ;
        }
      }

      LinAlg_AddVectorVector(&DofData_P->CurrentSolution->x, &DofData_P->dx, 
			     &DofData_P->CurrentSolution->x) ;
      
      Flag_CPU = 1 ;
      break ;

      /*  -->  S o l v e J a c _ A d a p t R e l a x  */
      /*  ------------------------------------------  */

    case OPERATION_SOLVEJACADAPTRELAX :
      /*  get increment dx by solving : J(xn) dx = b(xn) - A(xn) xn */
      Flag_Jac = 1 ;
      Init_OperationOnSystem("SolveJacAdaptRelax",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;

      if(DofData_P->Flag_Init[0] < 2)
	Msg(GERROR, "Jacobian system not initialized (missing GenerateJac?)");

      LinAlg_AddMatrixMatrix(&DofData_P->Jac, &DofData_P->A, &DofData_P->Jac) ;
      LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, &DofData_P->res) ;
      LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res) ;

      LinAlg_DummyVector(&DofData_P->res) ;

      LinAlg_Solve(&DofData_P->Jac, &DofData_P->res, &DofData_P->Solver, &DofData_P->dx) ;

      Msg(RESOURCES, "");
      
      /* save CurrentSolution */
      LinAlg_CreateVector(&x_Save, &DofData_P->Solver, DofData_P->NbrDof,
			  DofData_P->NbrPart, DofData_P->Part) ;
      LinAlg_CopyVector(&DofData_P->CurrentSolution->x, &x_Save);