PetscIJacobianIMEX.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__ "PetscIJacobianIMEX"

PetscErrorCode PetscIJacobianIMEX(TS ts,        
                                  PetscReal t,  
                                  Vec Y,        
                                  Vec Ydot,     
                                  PetscReal a,  
                                  Mat A,        
                                  Mat B,        
                                  void *ctxt     )
{
  PETScContext* context = (PETScContext*) ctxt;

  PetscFunctionBegin;
  for (int ns = 0; ns < context->nsims; ns++) {
    ((SimulationObject*)context->simobj)[ns].solver.count_IJacobian++;
  }
  context->shift = a;
  context->waqt  = t;
  /* Construct preconditioning matrix */
  if (context->flag_use_precon) PetscComputePreconMatIMEX(B,Y,context);

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "PetscJacobianFunctionIMEX_JFNK"

PetscErrorCode PetscJacobianFunctionIMEX_JFNK(Mat Jacobian, 
                                              Vec Y,
                                              Vec F  )
{
  PETScContext* context(nullptr);

  PetscFunctionBegin;

  MatShellGetContext(Jacobian,&context);
  SimulationObject* sim = (SimulationObject*) context->simobj;
  int nsims = context->nsims;

  double normY;
  VecNorm(Y,NORM_2,&normY);

  if (normY < 1e-16) {

    /* F = 0 */
    VecZeroEntries(F);
    /* [J]Y = aY - F(Y) */
    VecAXPBY(F,context->shift,0,Y);

  } else {
  
    double epsilon =  context->jfnk_eps / normY;
    double t = context->waqt; /* current stage/step time */

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

      HyPar* solver = &(sim[ns].solver);
      MPIVariables* mpi = &(sim[ns].mpi);
      solver->count_IJacFunction++;
  
      int size = solver->npoints_local_wghosts;
  
      double *u       = solver->u;
      double *uref    = solver->uref;
      double *rhsref  = solver->rhsref;
      double *rhs     = solver->rhs;
  
      /* copy solution from PETSc vector */
      TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
      _ArrayAYPX_(uref,epsilon,u,size*solver->nvars);
      /* 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 );
  
      /* Evaluate hyperbolic, parabolic and source terms  and the RHS for U+dU */
      _ArraySetValue_(rhs,size*solver->nvars,0.0);
      if ((!strcmp(solver->SplitHyperbolicFlux,"yes")) && solver->flag_fdf_specified) {
        if (context->flag_hyperbolic_f == _IMPLICIT_) {
          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 == _IMPLICIT_) {
          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 == _IMPLICIT_) {
          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 == _IMPLICIT_) {
          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 == _IMPLICIT_) {
          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 == _IMPLICIT_) {
        solver->ParabolicFunction (solver->par,u,solver,mpi,t);
        _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
      }
      if (context->flag_source == _IMPLICIT_) {
        solver->SourceFunction (solver->source,u,solver,mpi,t);
        _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
      }
  
      _ArrayAXPY_(rhsref,-1.0,rhs,size*solver->nvars);
      /* Transfer RHS to PETSc vector */
      TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
    }

    /* [J]Y = aY - F(Y) */
    VecAXPBY(F,context->shift,(-1.0/epsilon),Y);

  }

  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "PetscJacobianFunctionIMEX_Linear"

PetscErrorCode PetscJacobianFunctionIMEX_Linear(Mat Jacobian, 
                                                Vec Y,
                                                Vec F )
{
  PETScContext* context(nullptr);

  PetscFunctionBegin;

  MatShellGetContext(Jacobian,&context);
  SimulationObject* sim = (SimulationObject*) context->simobj;
  int nsims = context->nsims;

  double normY;
  VecNorm(Y,NORM_2,&normY);

  if (normY < 1e-16) {

    /* F = 0 */
    VecZeroEntries(F);
    /* [J]Y = aY - F(Y) */
    VecAXPBY(F,context->shift,0,Y);

  } else {
  
    double t = context->waqt; /* current stage/step time */

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

      HyPar* solver = &(sim[ns].solver);
      MPIVariables* mpi = &(sim[ns].mpi);
      solver->count_IJacFunction++;
  
      int size = solver->npoints_local_wghosts;
  
      double *u       = solver->u;
      double *uref    = solver->uref;
      double *rhsref  = solver->rhsref;
      double *rhs     = solver->rhs;
  
      /* copy solution from PETSc vector */
      TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
      _ArrayAYPX_(uref,1.0,u,size*solver->nvars);
      /* 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 );
  
      /* Evaluate hyperbolic, parabolic and source terms  and the RHS for U+dU */
      _ArraySetValue_(rhs,size*solver->nvars,0.0);
      if ((!strcmp(solver->SplitHyperbolicFlux,"yes")) && solver->flag_fdf_specified) {
        if (context->flag_hyperbolic_f == _IMPLICIT_) {
          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 == _IMPLICIT_) {
          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 == _IMPLICIT_) {
          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 == _IMPLICIT_) {
          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 == _IMPLICIT_) {
          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 == _IMPLICIT_) {
        solver->ParabolicFunction (solver->par,u,solver,mpi,t);
        _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
      }
      if (context->flag_source == _IMPLICIT_) {
        solver->SourceFunction (solver->source,u,solver,mpi,t);
        _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
      }
  
      _ArrayAXPY_(rhsref,-1.0,rhs,size*solver->nvars);
      /* Transfer RHS to PETSc vector */
      TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
    }

    /* [J]Y = aY - F(Y) */
    VecAXPBY(F,context->shift,-1.0,Y);

  }

  PetscFunctionReturn(0);
}

#endif