PetscIFunctionIMEX.cpp File Reference

Compute the implicitly-treated part of the right-hand-side for IMEX time integration. More...

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

Go to the source code of this file.

Macros
#define	__FUNCT__   "PetscIFunctionIMEX"

Functions
PetscErrorCode	PetscIFunctionIMEX (TS ts, PetscReal t, Vec Y, Vec Ydot, Vec F, void *ctxt)

Detailed Description

Compute the implicitly-treated part of the right-hand-side for IMEX time integration.

Author

Debojyoti Ghosh

Definition in file PetscIFunctionIMEX.cpp.

Macro Definition Documentation

#define __FUNCT__   "PetscIFunctionIMEX"

Definition at line 16 of file PetscIFunctionIMEX.cpp.

Function Documentation

PetscErrorCode PetscIFunctionIMEX	(	TS	ts,
		PetscReal	t,
		Vec	Y,
		Vec	Ydot,
		Vec	F,
		void *	ctxt
	)

Compute the implicitly-treated part of the right-hand-side for the implicit-explicit (IMEX) time integration of the governing equations: The ODE, obtained after discretizing the governing PDE in space, is expressed as follows (for the purpose of IMEX time integration):

\begin{eqnarray} \frac {d{\bf U}}{dt} &=& {\bf F}\left({\bf U}\right) + {\bf G}\left({\bf U}\right), \\ \Rightarrow \dot{\bf U} - {\bf G}\left({\bf U}\right) &=& {\bf F}\left({\bf U}\right), \end{eqnarray}

where \({\bf F}\) is non-stiff and integrated in time explicitly, and \({\bf G}\) is stiff and integrated in time implicitly, and \({\bf U}\) represents the entire solution vector (state vector).

Note that \({\bf G}\left({\bf U}\right)\) represents all terms that the user has indicated to be integrated in time implicitly (PETScContext::flag_hyperbolic_f, PETScContext::flag_hyperbolic_df, PETScContext::flag_hyperbolic, PETScContext::flag_parabolic, and PETScContext::flag_source).

This function computes the left-hand-side of the above equation:

\begin{equation} \mathcal{G}\left(\dot{\bf U},{\bf U},t\right) = \dot{\bf U} - {\bf G}\left({\bf U}\right) \end{equation}

given \(\dot{\bf U}\) and \({\bf U}\).

See Also

PetscRHSFunctionIMEX()

Notes:

• Y and Ydot in the code are \({\bf U}\) and \(\dot{\bf U}\), respectively. PETsc denotes the state vector with \({\bf Y}\) in its time integrators.
• It is assumed that the reader is familiar with PETSc's implementation of IMEX time integrators, for example, TSARKIMEX (https://petsc.org/release/docs/manualpages/TS/TSARKIMEX.html).
• See https://petsc.org/release/docs/manualpages/TS/index.html for documentation on PETSc's time integrators.
• All functions and variables whose names start with Vec, Mat, PC, KSP, SNES, and TS are defined by PETSc. Refer to the PETSc documentation (https://petsc.org/release/docs/). Usually, googling with the function or variable name yields the specific doc page dealing with that function/variable.

Parameters

ts	The time integration object
t	Current solution time
Y	State vector (input)
Ydot	Time derivative of the state vector (input)
F	The computed function vector
ctxt	Object of type PETScContext

Definition at line 53 of file PetscIFunctionIMEX.cpp.

{
  PETScContext* context = (PETScContext*) ctxt;
  SimulationObject* sim = (SimulationObject*) context->simobj;
  int nsims = context->nsims;

  PetscFunctionBegin;

  context->waqt = t;

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

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

    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 == _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);
    }
  
    /* save a copy of the solution and RHS for use in IJacobian */
    _ArrayCopy1D_(u  ,solver->uref  ,(size*solver->nvars));
    _ArrayCopy1D_(rhs,solver->rhsref,(size*solver->nvars));
  
    /* Transfer RHS to PETSc vector */
    TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);

  }

  /* LHS = Ydot - F(u) */
  VecAYPX(F,-1.0,Ydot);

  PetscFunctionReturn(0);
}