solvingoperations.c

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

	  Operation_P->Case.Generate_MH_Moving_S.Period ;
	Current.DTime = 1./(double)Operation_P->Case.Generate_MH_Moving_S.NbrStep * 
	  Operation_P->Case.Generate_MH_Moving_S.Period ;
	Current.TimeStep = iTime;

	if (!iTime) {
	  Msg(INFO, "Generate_MH_Moving_Separate : probing for any degrees of freedom"); 
	  DofTree_MH_moving = Tree_Create(sizeof(struct Dof), fcmp_Dof) ;	

	  /* probing assembly */
	  MH_Moving_Matrix_probe = 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) ;
	  }
	  MH_Moving_Matrix_probe = 0;

	  DofList_MH_moving = Tree2List(DofTree_MH_moving) ;	
	  Tree_Delete(DofTree_MH_moving) ;	

	  NbrDof_MH_moving = List_Nbr(DofList_MH_moving) ;
	  Msg(INFO,"Generate_MH_Moving :  NbrDof = %d", NbrDof_MH_moving);

	  Dof_MH_moving = (struct Dof **)Malloc(NbrDof_MH_moving * sizeof(struct Dof *)) ;
	  NumDof_MH_moving = (int *)Malloc(NbrDof_MH_moving * sizeof(int)) ;

	  for (i = 0 ; i < NbrDof_MH_moving ; i++) {	  
	    Dof_P = (struct Dof*)List_Pointer(DofList_MH_moving,i) ;
	    if (Dof_P->Type != DOF_UNKNOWN) Msg(GERROR,"Dof_MH_moving not of type unknown !?");
	    NumDof_MH_moving[i] =  Dof_P->Case.Unknown.NumDof; 

	    if(!(Dof_MH_moving[i] = (struct Dof *)List_PQuery(Current.DofData->DofList, 
							      Dof_P, fcmp_Dof)))
	      Msg(GERROR, "Troubles") ;
	    for (k = 0 ; k < Current.NbrHar ; k++) { 
	      (Dof_MH_moving[i]+k)->Case.Unknown.NumDof = i*Current.NbrHar+k+1 ;
	    }	  
	  } /* if (!iTime) */

	  Msg(RESOURCES, "");
	  LinAlg_CreateSolver(&DofData_P->Solver_MH_moving, "MH_moving.par") ;
	  LinAlg_CreateMatrix(&DofData_P->A_MH_moving, &DofData_P->Solver_MH_moving,
			      NbrDof_MH_moving*Current.NbrHar, NbrDof_MH_moving*Current.NbrHar,
			      DofData_P->NbrPart, DofData_P->Part, DofData_P->Nnz) ;
	  LinAlg_CreateVector(&DofData_P->b_MH_moving, &DofData_P->Solver_MH_moving, 
			      NbrDof_MH_moving*Current.NbrHar,
			      DofData_P->NbrPart, DofData_P->Part) ;
	  LinAlg_ZeroMatrix(&DofData_P->A_MH_moving) ;
	  LinAlg_ZeroVector(&DofData_P->b_MH_moving) ;      
	}

	Msg(INFO, "Generate_MH_Moving_Separate : Step %d/%d (Time = %e  DTime %e)", 
	    (int)(Current.TimeStep+1), Operation_P->Case.Generate_MH_Moving_S.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. ;
	
	/* separate assembly */
	MH_Moving_Matrix_separate = 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) ;
	}
	MH_Moving_Matrix_separate = 0 ;
	
	
      } /*  for iTime */

      Msg(RESOURCES, "Full matrix assembly done");

      for (k = 0 ; k < Current.NbrHar ; k++) Free (MH_Moving_Matrix[k]) ;
      Free (MH_Moving_Matrix) ;
      MH_Moving_Matrix = NULL ; 
      
      Generate_Group = NULL;
      
      for (i = 0 ; i < NbrDof_MH_moving ; i++) {	  
	for (k = 0 ; k < Current.NbrHar ; k++)	
	  (Dof_MH_moving[i]+k)->Case.Unknown.NumDof = NumDof_MH_moving[i] + k ;
      }
      
      LinAlg_CreateMatrix(&DofData_P->A_MH_moving2, &DofData_P->Solver,
			  DofData_P->NbrDof, DofData_P->NbrDof,
			  DofData_P->NbrPart, DofData_P->Part, DofData_P->Nnz) ;
      LinAlg_CreateVector(&DofData_P->b_MH_moving2, &DofData_P->Solver, Current.DofData->NbrDof,
			  DofData_P->NbrPart, DofData_P->Part) ;
      LinAlg_ZeroMatrix(&DofData_P->A_MH_moving2) ;
      LinAlg_ZeroVector(&DofData_P->b_MH_moving2) ;      

      Msg(RESOURCES, "");
      
      nnz__=0;
      for (i = 0 ; i < NbrDof_MH_moving ; i++) {
	for (k = 0 ; k < Current.NbrHar ; k++) {	  
	  row_old = Current.NbrHar*i+k ;
	  row_new = NumDof_MH_moving[i]+k-1 ;
	  LinAlg_GetDoubleInVector(&d, &DofData_P->b_MH_moving,  row_old) ;
	  LinAlg_SetDoubleInVector( d, &DofData_P->b_MH_moving2, row_new) ;
	  for (j = 0 ; j < NbrDof_MH_moving ; j++) {
	    for (l = 0 ; l < Current.NbrHar ; l++) {	  
	      col_old = Current.NbrHar*j+l ;
	      col_new = NumDof_MH_moving[j]+l-1 ;
	      
	      /* LinAlg_GetDoubleInMatrix(&d, &DofData_P->A_MH_moving, i, j) ; */
#if defined(HAVE_SPARSKIT)
	      d = DofData_P->A_MH_moving.M.F.a[NbrDof_MH_moving*Current.NbrHar*col_old+row_old]; 
	      aii = DofData_P->A_MH_moving.M.F.a[NbrDof_MH_moving*Current.NbrHar*row_old+row_old]; 
	      ajj = DofData_P->A_MH_moving.M.F.a[NbrDof_MH_moving*Current.NbrHar*col_old+col_old]; 
#else
	      aii = ajj = 0.;
	      Msg(GERROR, "FIXME: Generate_MH_Moving works only with Sparskit");
#endif
	      if(d*d > 1e-12 * aii*ajj  && 
		 ( (DummyDof[row_new]==0 && DummyDof[col_new] == 0) || (row_new == col_new) ) ){ 
		LinAlg_AddDoubleInMatrix(d, &DofData_P->A_MH_moving2, col_new, row_new) ;
		nnz__++;
	      }
	    }
	  }
	}
      }
      printf("Matrix converted : nnz in MH_moving %d \n", nnz__);
#if defined(HAVE_SPARSKIT)
      Free(DofData_P->A_MH_moving.M.F.a);
#endif
      Current.DTime = 0.;
      Msg(RESOURCES, "");
      DofData_P->DummyDof = DummyDof ;
      break;
      
    case OPERATION_DUMMYDOFS :    
      Init_OperationOnSystem("DummyDofs",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
                             &DefineSystem_P, &DofData_P, Resolution2_P) ;
      Msg(RESOURCES, "");
      Dof_GetDummies(DefineSystem_P, DofData_P);
      Msg(RESOURCES, "");
      break ;

    case OPERATION_ADD_MH_MOVING :
      LinAlg_AddMatrixMatrix(&DofData_P->A, &DofData_P->A_MH_moving2,&DofData_P->A) ;
      /* LinAlg_AddVectorVector(&DofData_P->b, &DofData_P->b_MH_moving2,&DofData_P->b) ; */
      Msg(RESOURCES, "");
      break ;

      /*  -->  S a v e S o l u t i o n E x t e n d e d M H             */
      /*  -----------------------------------------------------------  */

    case OPERATION_SAVESOLUTIONEXTENDEDMH :
      if (Current.NbrHar == 1) { 
	Msg(WARNING, "ExtendSolutionMH can only to be used with multi-harmonics") ;
	break ;
      }
      else if (!List_Nbr(DofData_P->Solutions)) { 
	Msg(WARNING, "No solution available for ExtendSolutionMH");
	break ;
      }
      else if (List_Nbr(DofData_P->Solutions) > 1) {
	Msg(WARNING, "Only last solution will be extended mult-harmonically and saved");
      }

      Init_OperationOnSystem("SaveSolutionExtendedMH",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
			     &DefineSystem_P, &DofData_P, Resolution2_P) ;
      strcpy(FileName_exMH, Name_Generic) ;
      strcat(FileName_exMH, Operation_P->Case.SaveSolutionExtendedMH.ResFile) ;
      strcat(FileName_exMH, ".res") ;
      Dof_WriteFileRES0(FileName_exMH, Flag_BIN) ;      
      Dof_WriteFileRES_ExtendMH(FileName_exMH, DofData_P, Flag_BIN,  
				Current.NbrHar + 
				2*Operation_P->Case.SaveSolutionExtendedMH.NbrFreq);

      Msg(DIRECT, "          > '%s'  (%d to %d frequencies)", FileName_exMH, 
	  Current.NbrHar/2, Current.NbrHar/2 + Operation_P->Case.SaveSolutionExtendedMH.NbrFreq) ;

      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 M H T o T i m e                 */
      /*  -----------------------------------------------------------  */

    case OPERATION_SAVESOLUTIONMHTOTIME :
      if (Current.NbrHar == 1) { 
	Msg(WARNING, "SaveSolutionMHtoTime can only to be used with multi-harmonics") ;
	break ;
      } 
      else if (!List_Nbr(DofData_P->Solutions)) { 
	Msg(WARNING, "No solution available for SaveSolutionMHtoTime");
	break ;
      }
      else if (List_Nbr(DofData_P->Solutions) > 1) {
	Msg(WARNING, "Only last mult-harmonic solution will be saved for time X");
      }

      Init_OperationOnSystem("SaveSolutionMHtoTime",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
			     &DefineSystem_P, &DofData_P, Resolution2_P) ;
      strcpy(FileName_exMH, Name_Generic) ;
      strcat(FileName_exMH, Operation_P->Case.SaveSolutionMHtoTime.ResFile) ;
      strcat(FileName_exMH, ".res") ;
      Dof_WriteFileRES0(FileName_exMH, Flag_BIN) ;      
      Dof_WriteFileRES_MHtoTime(FileName_exMH, DofData_P, Flag_BIN, 
				Operation_P->Case.SaveSolutionMHtoTime.Time) ;

      Msg(DIRECT, "      > '%s'  (time = %e)", FileName_exMH, 
	  Operation_P->Case.SaveSolutionMHtoTime.Time) ;

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

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

    case OPERATION_READSOLUTION :
      Init_OperationOnSystem("ReadSolution",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
			     &DefineSystem_P, &DofData_P, Resolution2_P) ;
      i = 0 ;
      while(Name_ResFile[i]){
	Msg(LOADING, "Processing data '%s'", Name_ResFile[i]) ;
	Dof_OpenFile(DOF_TMP, Name_ResFile[i], "rb");
	Dof_ReadFileRES(NULL, DofData_P, DofData_P->Num, &Current.Time, &Current.TimeImag,
			&Current.TimeStep) ;
	Dof_CloseFile(DOF_TMP);
	i++ ;
      }
      if(!List_Nbr(DofData_P->Solutions))
	Msg(GERROR, "No valid data found for ReadSolution[%s]", DefineSystem_P->Name);
	
      DofData_P->CurrentSolution = (struct Solution*)
	List_Pointer(DofData_P->Solutions, List_Nbr(DofData_P->Solutions)-1) ;	
      DofData_P->CurrentSolution->TimeFunctionValues = Get_TimeFunctionValues(DofData_P) ;
      break ;

      /*  -->  S a v e M e s h                        */
      /*  ------------------------------------------  */

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

      if(Operation_P->Case.SaveMesh.FileName[0] == '/' || 
	 Operation_P->Case.SaveMesh.FileName[0] == '\\'){
	strcpy(FileName,Operation_P->Case.SaveMesh.FileName);
      } else {
	strcpy(FileName, Name_Path);
	strcat(FileName, Operation_P->Case.SaveMesh.FileName);
      }

      if (Operation_P->Case.SaveMesh.ExprIndex >= 0) {
	Get_ValueOfExpressionByIndex(Operation_P->Case.SaveMesh.ExprIndex,
				     NULL, 0., 0., 0., &Value) ;
	sprintf(FileName, FileName, Value.Val[0]);
      } 

      Geo_SaveMesh(Current.GeoData, 
		   ((struct Group*)
		    List_Pointer(Problem_S.Group, 
				 Operation_P->Case.SaveMesh.GroupIndex))->InitialList,
		   FileName) ;
      break ;

      /*  -->  T r a n s f e r S o l u t i o n        */
      /*  ------------------------------------------  */

    case OPERATION_TRANSFERSOLUTION :
      Init_OperationOnSystem("TransferSolution",
			     Resolution_P, Operation_P, DofData_P0, GeoData_P0,
			     &DefineSystem_P, &DofData_P, Resolution2_P) ;
      
      if(Resolution2_P){ /* pre-resolution */
	DofData2_P = DofData2_P0 + DefineSystem_P->DestinationSystemIndex ;
	Dof_TransferDof(DofData_P, &DofData2_P);
      }
      else{ 
	/* a changer!!! Il faut se mettre d'accord sur ce que doit faire 
	   Dof_TransferDof. Ceci sert a transferer la derniere solution d'un 
	   DofData dans un autre (ds la meme resolution), base sur le meme 
	   espace fonctionnel. */
	DofData2_P = DofData_P0 + DefineSystem_P->DestinationSystemIndex ;
	
	if(DofData_P->NbrAnyDof != DofData2_P->NbrAnyDof)
	  Msg(GERROR, "Dimensions do not match for TransferSolution");

	Solution_S.TimeStep = (int)Current.TimeStep ;
	Solution_S.Time = Current.Time ;
	Solution_S.TimeImag = Current.TimeImag ;
	Solution_S.TimeFunctionValues = Get_TimeFunctionValues(DofData2_P) ;
	Solution_S.SolutionExist = 1 ;
        LinAlg_CreateVector(&Solution_S.x, &DofData2_P->Solver, DofData2_P->NbrDof,
			    DofData2_P->NbrPart, DofData2_P->Part) ;
	LinAlg_ZeroVector(&Solution_S.x) ;

	if (List_Nbr(DofData_P->Solutions)) {
	  
	  Solution_P = (struct Solution *)List_Pointer(DofData_P->Solutions,
						       List_Nbr(DofData_P->Solutions)-1) ;
	  for(i=0 ; i<DofData_P->NbrAnyDof ; i++){
	    Dof = *(struct Dof *)List_Pointer(DofData_P->DofList, i) ;
	    if(Dof.Type == DOF_UNKNOWN){
	      LinAlg_GetScalarInVector(&tmp, &Solution_P->x, Dof.Case.Unknown.NumDof-1) ;

	      if((Dof_P = (struct Dof*)List_PQuery(DofData2_P->DofList, &Dof, fcmp_Dof))){
		LinAlg_SetScalarInVector(&tmp, &Solution_S.x, Dof_P->Case.Unknown.NumDof-1) ;
		Dof_P->Type = DOF_UNKNOWN ;
	      }
	      else{
		Msg(WARNING, "Unknown Dof in TransferSolution") ;
	      }
	    }
	    else{
	      Msg(WARNING, "Trying to transfer a non symmetrical Dof");
	    }
	  }

	  if (!DofData2_P->Solutions)
	    DofData2_P->Solutions = List_Create( 20, 20, sizeof(struct Solution)) ;

	  List_Add(DofData2_P->Solutions, &Solution_S) ;      
	  DofData2_P->CurrentSolution = (struct Solution*)
	    List_Pointer(DofData2_P->Solutions, List_Nbr(DofData2_P->Solutions)-1) ;
	}
      }
      break ;

      /*  -->  E v a l u a t e                        */
      /*  ------------------------------------------  */

    case OPERATION_EVALUATE :
      Get_ValueOfExpressionByIndex(Operation_P->Case.Evaluate.ExpressionIndex,
				   NULL, 0., 0., 0., &Value) ;
      break ;