PetscRHSFunctionIMEX.cpp

Go to the documentation of this file.

#ifdef with_petsc

#include <stdlib.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mpivars_cpp.h>
#include <simulation_object.h>
#include <petscinterface.h>

#undef __FUNCT__
#define __FUNCT__ "PetscRHSFunctionIMEX"

PetscErrorCode PetscRHSFunctionIMEX(  TS        ts,   
                                      PetscReal t,    
                                      Vec       Y,    
                                      Vec       F,    
                                      void      *ctxt  )
{
  PETScContext* context = (PETScContext*) ctxt;
  SimulationObject* sim = (SimulationObject*) context->simobj;
  int nsims = context->nsims;

  context->waqt = t;

  PetscFunctionBegin;
  for (int ns = 0; ns < nsims; ns++) {

    HyPar* solver = &(sim[ns].solver);
    MPIVariables* mpi = &(sim[ns].mpi);
  
    solver->count_RHSFunction++;

    int size = solver->npoints_local_wghosts;
    double *u = solver->u;
    double *rhs = solver->rhs;
  
    /* copy solution from PETSc vector */
    TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
    /* apply boundary conditions and exchange data over MPI interfaces */
    solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
    MPIExchangeBoundariesnD(solver->ndims,solver->nvars,solver->dim_local,
                                 solver->ghosts,mpi,u);
  
    /* initialize right-hand side to zero */
    _ArraySetValue_(rhs,size*solver->nvars,0.0);
  
    /* Evaluate hyperbolic, parabolic and source terms  and the RHS */
    if ((!strcmp(solver->SplitHyperbolicFlux,"yes")) && solver->flag_fdf_specified) {
      if (context->flag_hyperbolic_f == _EXPLICIT_) {
        solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FdFFunction,solver->UpwindFdF);  
        _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
      } 
      if (context->flag_hyperbolic_df == _EXPLICIT_) {
        solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF);
        _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
      }
    } else if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
      if (context->flag_hyperbolic_f == _EXPLICIT_) {
        solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FFunction,solver->Upwind);  
        _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
        solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF); 
        _ArrayAXPY_(solver->hyp, 1.0,rhs,size*solver->nvars);
      } 
      if (context->flag_hyperbolic_df == _EXPLICIT_) {
        solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF); 
        _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
      }
    } else {
      if (context->flag_hyperbolic == _EXPLICIT_) {
        solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FFunction,solver->Upwind);  
        _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
      }
    }
    if (context->flag_parabolic == _EXPLICIT_) {
      solver->ParabolicFunction (solver->par,u,solver,mpi,t);                        
      _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
    }
    if (context->flag_source == _EXPLICIT_) {
      solver->SourceFunction    (solver->source,u,solver,mpi,t);                     
      _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
    }
  
    /* Transfer RHS to PETSc vector */
    TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
  }

  PetscFunctionReturn(0);
}

#endif