solvingoperations.c

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

      Flag_RHS = 1;  
      /* MHJacNL-terms don't contribute to the RHS and residu, and are thus disregarded */

      Error_Prev = 1e99 ;  Frelax_Opt = 1. ;

      if (!(NbrSteps_relax = List_Nbr(Operation_P->Case.SolveJac_AdaptRelax.Factor_L)))
	  Msg(GERROR, "No factors provided for Adaptive Relaxation");

      for( istep = 0 ; istep < NbrSteps_relax ; istep++ ){  
		
	List_Read(Operation_P->Case.SolveJac_AdaptRelax.Factor_L, istep, &Frelax);

	/* new trial solution = x + Frelax * dx */
	LinAlg_CopyVector(&x_Save, &DofData_P->CurrentSolution->x);
	LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x, &DofData_P->dx, 
					 Frelax, &DofData_P->CurrentSolution->x);
	/* LinAlg_PrintVector(stdout, &DofData_P->CurrentSolution->x); */

	/* calculate residual with trial solution */
	ReGenerate_System(DefineSystem_P, DofData_P, DofData_P0) ;
	LinAlg_ProdMatrixVector(&DofData_P->A, &DofData_P->CurrentSolution->x, &DofData_P->res) ;
	LinAlg_SubVectorVector(&DofData_P->b, &DofData_P->res, &DofData_P->res) ;

	/* check whether norm of residual is smaller than previous ones */
	LinAlg_VectorNorm2(&DofData_P->res, &Norm);
	LinAlg_GetVectorSize(&DofData_P->res, &N);
	Norm /= (double)N;
	Msg(INFO, " adaptive relaxation : factor = %8f   Norm residual = %10.4e", Frelax, Norm) ;

	if (Norm < Error_Prev) {
	  Error_Prev = Norm;
	  Frelax_Opt = Frelax;
	} else if ( !Operation_P->Case.SolveJac_AdaptRelax.CheckAll && istep > 0 ) break ;

      }

      Msg(INFO, " => optimal relaxation factor = %f", Frelax_Opt) ;

      /*  solution = x + Frelax_Opt * dx */
      LinAlg_CopyVector(&x_Save, &DofData_P->CurrentSolution->x);
      LinAlg_AddVectorProdVectorDouble(&DofData_P->CurrentSolution->x, &DofData_P->dx, 
				       Frelax_Opt, &DofData_P->CurrentSolution->x);

      MeanError = Error_Prev ; 
      Msg(BIGINFO, "Mean error: %.3e  (after %d iteration%s)", 
	  MeanError, (int)Current.Iteration, ((int)Current.Iteration==1)?"":"s") ;

      Current.RelativeDifference = MeanError;
      Flag_CPU = 1 ;
      Flag_RHS = 0 ;
      LinAlg_DestroyVector(&x_Save);
      break ;

      /*  -->  Lanczos                                */
      /*  ------------------------------------------  */

    case OPERATION_LANCZOS :
      Init_OperationOnSystem("Lanczos",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;
      Lanczos(DofData_P, Operation_P->Case.Lanczos.Size, 
	      Operation_P->Case.Lanczos.Save, Operation_P->Case.Lanczos.Shift) ;
      Flag_CPU = 1 ;
      break ;

      /*  -->  EigenSolve                             */
      /*  ------------------------------------------  */

    case OPERATION_EIGENSOLVE :
      Init_OperationOnSystem("EigenSolve",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;
      EigenSolve(DofData_P, Operation_P->Case.EigenSolve.NumEigenvalues, 
		 Operation_P->Case.EigenSolve.Shift_r,
		 Operation_P->Case.EigenSolve.Shift_i) ;
      Flag_CPU = 1 ;
      break ;

      /*  -->  EigenSolveJac                             */
      /*  ------------------------------------------  */

    case OPERATION_EIGENSOLVEJAC :
      Init_OperationOnSystem("EigenSolveJac",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;
      EigenSolve(DofData_P, Operation_P->Case.EigenSolve.NumEigenvalues, 
		 Operation_P->Case.EigenSolve.Shift_r,
		 Operation_P->Case.EigenSolve.Shift_i) ;
      /* Insert intelligent convergence test here :-) */
      Current.RelativeDifference = 1.0 ;
      Flag_CPU = 1 ;
      break ;

      /*  -->  Perturbation                           */
      /*  ------------------------------------------  */

    case OPERATION_PERTURBATION :
      Init_OperationOnSystem("Perturbation",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;
      /*
      Perturbation(DofData_P,
		   DofData_P0+Operation_P->Case.Perturbation.DefineSystemIndex2,
		   DofData_P0+Operation_P->Case.Perturbation.DefineSystemIndex3,
		   Operation_P->Case.Perturbation.Size, 
	           Operation_P->Case.Perturbation.Save,
		   Operation_P->Case.Perturbation.Shift,
		   Operation_P->Case.Perturbation.PertFreq) ;
      */
      Flag_CPU = 1 ;
      break ;

      /*  -->  I n i t S o l u t i o n                */
      /*  ------------------------------------------  */

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

      if(Flag_RESTART){
        if (!DofData_P->Solutions)
          Msg(GERROR, "No solution to restart the computation");

        for(i = 0 ; i < DofData_P->NbrAnyDof ; i++){
          Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, i) ;
          if(Dof_P->Type == DOF_UNKNOWN_INIT) Dof_P->Type = DOF_UNKNOWN ;
        }

	for(i = 0; i < List_Nbr(DofData_P->Solutions); i++){
	  Solution_P = (struct Solution*)List_Pointer(DofData_P->Solutions, i);
	  Solution_P->TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ;
	  /* The last solution is the current one */	  
	  if(i == List_Nbr(DofData_P->Solutions) - 1)
	    DofData_P->CurrentSolution = Solution_P;
	}
	RES0 = (int)Current.TimeStep ;
      }
      else{
        if (!DofData_P->Solutions)
          DofData_P->Solutions = List_Create( 20, 20, sizeof(struct Solution)) ;
        
	Solution_S.TimeStep = (int)Current.TimeStep ;
        Solution_S.Time = Current.Time ;
        Solution_S.TimeImag = Current.TimeImag ;
        Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ;

	Solution_S.SolutionExist = 1 ;
        LinAlg_CreateVector(&Solution_S.x, &DofData_P->Solver, DofData_P->NbrDof,
			    DofData_P->NbrPart, DofData_P->Part) ;

	/* The last solution, if any, initializes the current one.
	   Otherwise a null solution is used. 
	   a revoir qd les conditions initiales multiples seront mieux traitees */
	if (List_Nbr(DofData_P->Solutions)) {
	  LinAlg_CopyVector(&((struct Solution *)
			      List_Pointer(DofData_P->Solutions,
					   List_Nbr(DofData_P->Solutions)-1))->x,
			    &Solution_S.x) ;
	}
	else {
	  LinAlg_ZeroVector(&Solution_S.x) ;
	}

	for(i=0 ; i<DofData_P->NbrAnyDof ; i++){
	  Dof_P = (struct Dof *)List_Pointer(DofData_P->DofList, i) ;
	  if(Dof_P->Type == DOF_UNKNOWN_INIT){ /* Init values loaded */
	    Dof_P->Type = DOF_UNKNOWN ;
	    LinAlg_SetScalarInVector
	      (&Dof_P->Val, &Solution_S.x, Dof_P->Case.Unknown.NumDof-1) ;
	  }
	}
	
	List_Add(DofData_P->Solutions, &Solution_S) ;

	DofData_P->CurrentSolution = (struct Solution*)
	  List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ;
      }
      break ;

      /*  -->  S a v e S o l u t i o n                */
      /*  ------------------------------------------  */

    case OPERATION_SAVESOLUTION :
      Init_OperationOnSystem("SaveSolution",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
			     &DefineSystem_P, &DofData_P, Resolution2_P) ;
      strcpy(ResName, Name_Generic) ;
      if(!Flag_SPLIT){
	strcat(ResName, ".res") ;
	if(RES0 < 0){
	  Dof_WriteFileRES0(ResName, Flag_BIN) ;
	  RES0 = 1 ;
	}
      }
      else{
	strcat(ResName, "-") ;
	sprintf(ResNum, "%d.res", (int)Current.TimeStep) ;
	for(i = 0 ; i < 5+4-(int)strlen(ResNum) ; i++) strcat(ResName, "0") ;
	strcat(ResName, ResNum) ;
	if(RES0 != (int)Current.TimeStep){
	  Dof_WriteFileRES0(ResName, Flag_BIN) ;
	  RES0 = (int)Current.TimeStep ;
	}
      }
      Dof_WriteFileRES(ResName, DofData_P, Flag_BIN, Current.Time, Current.TimeImag,
		       (int)Current.TimeStep) ;
      break ;

      /*  -->  S a v e S o l u t i o n s              */
      /*  ------------------------------------------  */

    case OPERATION_SAVESOLUTIONS :
      Init_OperationOnSystem("SaveSolutions",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
			     &DefineSystem_P, &DofData_P, Resolution2_P) ;
      strcpy(ResName, Name_Generic) ;
      strcat(ResName, ".res") ;
      if(RES0 < 0){	
	Dof_WriteFileRES0(ResName, Flag_BIN) ;
	RES0 = 1 ;
      }
      for(i=0 ; i<List_Nbr(DofData_P->Solutions) ; i++){
	DofData_P->CurrentSolution = (struct Solution*)
	  List_Pointer(DofData_P->Solutions, i) ;
	if (!DofData_P->CurrentSolution->SolutionExist)
	  Msg(WARNING, "Solution #%d doesn't exist anymore: skipping", i) ;
	else
	  Dof_WriteFileRES(ResName, DofData_P, Flag_BIN, 
			   DofData_P->CurrentSolution->Time, 
			   DofData_P->CurrentSolution->TimeImag, i) ;
      }
      break ;

      /*  -->  M o v i n g   B a n d                  */
      /*  ------------------------------------------  */

    case OPERATION_INIT_MOVINGBAND2D :
      Init_MovingBand2D( (struct Group *)
			 List_Pointer(Problem_S.Group, 
				      Operation_P->Case.Init_MovingBand2D.GroupIndex)) ;
      break ;

    case OPERATION_MESH_MOVINGBAND2D :
      Mesh_MovingBand2D( (struct Group *)
			 List_Pointer(Problem_S.Group, 
				      Operation_P->Case.Mesh_MovingBand2D.GroupIndex)) ;
      break ;


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

      if(gSCALAR_SIZE == 2)
	Msg(GERROR, "FIXME: Generate_MH_Moving will not work in complex arithmetic");
      
      Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ;

      Generate_Group = (struct Group *) 
	List_Pointer(Problem_S.Group, Operation_P->Case.Generate_MH_Moving.GroupIndex) ;

      MH_Moving_Matrix = (double **) Malloc(Current.NbrHar*sizeof(double *)) ;
      for (k = 0 ; k < Current.NbrHar ; k++)
	MH_Moving_Matrix[k] = (double *) Malloc(Current.NbrHar*sizeof(double)) ;

      if (! (Val_Pulsation = Current.DofData->Val_Pulsation))
	Msg(GERROR, "Generate_MH_moving can only be used for harmonic problems");

      for (k = 0 ; k < Current.NbrHar ; k++)
	for (l = 0 ; l < Current.NbrHar ; l++) 
	  hop[k][l] = 0.;

      MH_Moving_Matrix_simple = 1;
      
      for (iTime = 0 ; iTime < Operation_P->Case.Generate_MH_Moving.NbrStep ; iTime++) {
      
	Current.Time = (double)iTime/(double)Operation_P->Case.Generate_MH_Moving.NbrStep * 
	  Operation_P->Case.Generate_MH_Moving.Period ;
	Current.DTime = 1./(double)Operation_P->Case.Generate_MH_Moving.NbrStep * 
	  Operation_P->Case.Generate_MH_Moving.Period ;
	Current.TimeStep = iTime;

	Msg(INFO, "Generate_MH_Moving : Step %d/%d (Time = %e  DTime %e)", 
	    (int)(Current.TimeStep+1), Operation_P->Case.Generate_MH_Moving.NbrStep, 
	    Current.Time, Current.DTime) ;

	Treatment_Operation(Resolution_P, Operation_P->Case.Generate_MH_Moving.Operation, 
			    DofData_P0, GeoData_P0, NULL, NULL) ;

	for (k = 0 ; k < Current.NbrHar ; k++)
	  for (l = 0 ; l < Current.NbrHar ; l++) {
	    if (Val_Pulsation[k/2]) DCfactor = 2. ; else DCfactor = 1. ; 
	    MH_Moving_Matrix[k][l] = DCfactor / 
	      (double)Operation_P->Case.Generate_MH_Moving.NbrStep *
	      ( fmod(k,2) ? -sin(Val_Pulsation[k/2]*Current.Time) :
		cos(Val_Pulsation[k/2]*Current.Time) ) *
	      ( fmod(l,2) ? -sin(Val_Pulsation[l/2]*Current.Time) : 
		cos(Val_Pulsation[l/2]*Current.Time) ) ;
	    hop[k][l] += MH_Moving_Matrix[k][l] ;
	  }

	for (k = 0 ; k < Current.NbrHar/2 ; k++)
	  if (!Val_Pulsation[k]) MH_Moving_Matrix[2*k+1][2*k+1] = 1. ;
	
	for (i = 0 ; i < Nbr_Formulation ; i++) {
	  List_Read(DefineSystem_P->FormulationIndex, i, &Index_Formulation) ;
	  Formulation_P = (struct Formulation*)
	    List_Pointer(Problem_S.Formulation, Index_Formulation) ;
	  Treatment_Formulation(Formulation_P) ;
	}
	
      }

      Current.TimeStep = 0;
      Current.Time = 0.;

      for (k = 0 ; k < Current.NbrHar ; k++) Free (MH_Moving_Matrix[k]) ;
      Free (MH_Moving_Matrix) ;
      MH_Moving_Matrix = NULL ; 
      MH_Moving_Matrix_simple = 0 ;

      Generate_Group = NULL;

      Msg(RESOURCES, "");
      break ;

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

      Nbr_Formulation = List_Nbr(DefineSystem_P->FormulationIndex) ;

      Generate_Group = (struct Group *)
	List_Pointer(Problem_S.Group, Operation_P->Case.Generate_MH_Moving_S.GroupIndex) ;

      MH_Moving_Matrix = (double **) Malloc(Current.NbrHar*sizeof(double *)) ;
      for (k = 0 ; k < Current.NbrHar ; k++)
	MH_Moving_Matrix[k] = (double *) Malloc(Current.NbrHar*sizeof(double)) ;

      if (! (Val_Pulsation = Current.DofData->Val_Pulsation))
	Msg(GERROR, "Generate_MH_moving can only be used for harmonic problems");

      for (k = 0 ; k < Current.NbrHar ; k++)
	for (l = 0 ; l < Current.NbrHar ; l++) 
	  hop[k][l] = 0.;

      DummyDof = DofData_P->DummyDof ;
      DofData_P->DummyDof = NULL ;
      
      for (iTime = 0 ; iTime < Operation_P->Case.Generate_MH_Moving_S.NbrStep ; iTime++) {
      
	Current.Time = (double)iTime/(double)Operation_P->Case.Generate_MH_Moving_S.NbrStep *