solutionsequence_fct.cpp File Reference

: numerical core of flux corrected transport solution More...

#include "../toolkits/toolkits.h"
#include "../classes/classes.h"
#include "../shared/shared.h"
#include "../modules/modules.h"
#include "../analyses/analyses.h"

Go to the source code of this file.

Macros
#define	USEPENALTYMETHOD   false

Functions
void	solutionsequence_fct (FemModel *femmodel)

Detailed Description

: numerical core of flux corrected transport solution

Definition in file solutionsequence_fct.cpp.

Macro Definition Documentation

USEPENALTYMETHOD

#define USEPENALTYMETHOD   false

Definition at line 10 of file solutionsequence_fct.cpp.

Function Documentation

solutionsequence_fct()

void solutionsequence_fct

(

FemModel *

femmodel

)

Definition at line 388 of file solutionsequence_fct.cpp.

                                             {/*{{{*/
 
   /*intermediary: */
   IssmDouble           theta,deltat,dmax;
   int                  configuration_type,analysis_type;
   Matrix<IssmDouble>*  K  = NULL;
   Matrix<IssmDouble>*  Mc = NULL;
   Vector<IssmDouble>*  p  = NULL;
   Vector<IssmDouble>*  ug = NULL;
   Vector<IssmDouble>*  uf = NULL;
   MasstransportAnalysis* manalysis = NULL;
   DamageEvolutionAnalysis* danalysis = NULL;
 
   /*Recover parameters: */
   femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
   femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
   femmodel->UpdateConstraintsx();
   theta = 0.5; //0.5==Crank-Nicolson, 1==Backward Euler
 
   /*Create analysis*/
   femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
   switch(analysis_type){
      case MasstransportAnalysisEnum:
         manalysis = new MasstransportAnalysis();
         manalysis->MassMatrix(&Mc,femmodel);
         manalysis->FctKMatrix(&K,NULL,femmodel);
         manalysis->FctPVector(&p,femmodel);
         break;
      case DamageEvolutionAnalysisEnum:
         danalysis = new DamageEvolutionAnalysis();
         danalysis->MassMatrix(&Mc,femmodel);
         danalysis->FctKMatrix(&K,NULL,femmodel);
         break;
      default: _error_("analysis type " << EnumToStringx(analysis_type) << " not supported for FCT\n");
   }
   delete manalysis;
   delete danalysis;
 
   #ifdef _HAVE_PETSC_
 
   /*Convert matrices to PETSc matrices*/
   Mat D_petsc  = NULL;
   Mat LHS      = NULL;
   Vec RHS      = NULL;
   Vec u        = NULL;
   Vec udot     = NULL;
   Mat K_petsc  = K->pmatrix->matrix;
   Mat Mc_petsc = Mc->pmatrix->matrix;
   Vec p_petsc  = p->pvector->vector;
 
   /*Get previous solution u^n*/
   GetSolutionFromInputsx(&ug,femmodel);
   Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
   delete ug;
 
   /*Compute lumped mass matrix*/
   Vec Ml_petsc = NULL;
   VecDuplicate(uf->pvector->vector,&Ml_petsc);
   MatGetRowSum(Mc_petsc,Ml_petsc);
 
   /*Create D Matrix*/
   CreateDMatrix(&D_petsc,K_petsc);
 
   /*Create LHS: [ML - theta*detlat *(K+D)^n+1]*/
   CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
 
   /*Create RHS: [ML + (1 - theta) deltaT L^n] u^n */
   CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,p_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
   delete uf;
   delete p;
 
   /*Go solve lower order solution*/
   femmodel->profiler->Start(SOLVER);
   SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters); 
   femmodel->profiler->Stop(SOLVER);
   MatFree(&LHS);
   VecFree(&RHS);
 
   /*Richardson to calculate udot*/
   RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
   delete K;
 
   /*Serialize u and udot*/
   IssmDouble* udot_serial = NULL;
   IssmDouble* u_serial    = NULL;
   IssmDouble* ml_serial   = NULL;
   VecToMPISerial(&udot_serial,udot    ,IssmComm::GetComm());
   VecToMPISerial(&u_serial   ,u       ,IssmComm::GetComm());
   VecToMPISerial(&ml_serial  ,Ml_petsc,IssmComm::GetComm());
 
   /*Anti diffusive fluxes*/
   Vec         Fbar            = NULL;
   IssmDouble *Ri_minus_serial = NULL;
   IssmDouble *Ri_plus_serial  = NULL;
   IssmDouble *ulmin           = NULL;
   IssmDouble *ulmax           = NULL;
   femmodel->GetInputLocalMinMaxOnNodesx(&ulmin,&ulmax,u_serial);
   CreateRis(&Ri_plus_serial,&Ri_minus_serial,Mc_petsc,D_petsc,ml_serial,u,u_serial,udot_serial,ulmin,ulmax,deltat);
   CreateFbar(&Fbar,Ri_plus_serial,Ri_minus_serial,Mc_petsc,D_petsc,udot_serial,u_serial,u);
   xDelete<IssmDouble>(Ri_plus_serial);
   xDelete<IssmDouble>(Ri_minus_serial);
   xDelete<IssmDouble>(ulmin);
   xDelete<IssmDouble>(ulmax);
 
   /*Clean up*/
   MatFree(&D_petsc);
   delete Mc;
   xDelete<IssmDouble>(udot_serial);
   xDelete<IssmDouble>(u_serial);
   xDelete<IssmDouble>(ml_serial);
 
   /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
   UpdateSolution(u,udot,Ml_petsc,Fbar,deltat);
   uf =new Vector<IssmDouble>(u);
   VecFree(&u);
   VecFree(&Fbar);
   VecFree(&udot);
   VecFree(&Ml_petsc);
 
   /*Update Element inputs*/
   InputUpdateFromSolutionx(femmodel,uf); 
   delete uf;
 
   #else
   _error_("PETSc needs to be installed");
   #endif
}/*}}}*/