PetscIJacobianIMEX.cpp File Reference

Contains the functions required for Jacobian computations 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__   "PetscIJacobianIMEX"
#define	__FUNCT__   "PetscJacobianFunctionIMEX_JFNK"
#define	__FUNCT__   "PetscJacobianFunctionIMEX_Linear"

Functions
PetscErrorCode	PetscIJacobianIMEX (TS ts, PetscReal t, Vec Y, Vec Ydot, PetscReal a, Mat A, Mat B, void *ctxt)
PetscErrorCode	PetscJacobianFunctionIMEX_JFNK (Mat Jacobian, Vec Y, Vec F)
PetscErrorCode	PetscJacobianFunctionIMEX_Linear (Mat Jacobian, Vec Y, Vec F)

Detailed Description

Contains the functions required for Jacobian computations for IMEX time-integration.

Author

Debojyoti Ghosh

Definition in file PetscIJacobianIMEX.cpp.

Macro Definition Documentation

#define __FUNCT__   "PetscIJacobianIMEX"

Definition at line 237 of file PetscIJacobianIMEX.cpp.

#define __FUNCT__   "PetscJacobianFunctionIMEX_JFNK"

Definition at line 237 of file PetscIJacobianIMEX.cpp.

#define __FUNCT__   "PetscJacobianFunctionIMEX_Linear"

Definition at line 237 of file PetscIJacobianIMEX.cpp.

Function Documentation

PetscErrorCode PetscIJacobianIMEX	(	TS	ts,
		PetscReal	t,
		Vec	Y,
		Vec	Ydot,
		PetscReal	a,
		Mat	A,
		Mat	B,
		void *	ctxt
	)

Compute the Jacobian for implicit-explicit (IMEX) time integration of the governing equations: The ODE, obtained after discretizing the governing PDE in space, is expressed as follows:

\begin{equation} \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{equation}

where \({\bf F}\) is the spatially discretized non-stiff right-hand-side (treated explicitly), \({\bf G}\) is the spatially discretized stiff right-hand-side (treated implicitly), and \({\bf U}\) represents the entire solution vector. The Jacobian of the implicit part is thus given by:

\begin{equation} {\bf J} = \left[\alpha{\bf I} - \frac {\partial {\bf G}} {\partial {\bf U}} \right] \end{equation}

where \(\alpha\) is the shift coefficient (PETScContext::shift) of the time integration method.

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).

Matrix-free representation: The Jacobian is computed using a matrix-free approach, where the entire Jacobian is never assembled, computed, or stored. Instead, the action of the Jacobian on a vector is defined through functions (PetscJacobianFunctionIMEX_JFNK() for nonlinear systems, and PetscJacobianFunctionIMEX_Linear() for linear systems). This approach works well with iterative solvers. Thus, this function just does the following:

• Saves the PETScContext::shift ( \(\alpha\)) and PETScContext::waqt (current simulation time) to the application context (so that PetscJacobianFunctionIMEX_JFNK() or PetscJacobianFunctionIMEX_Linear() can access these values).
• If a preconditioner is being used, calls the function to compute the preconditioning matrix.

See Also

PetscIFunctionIMEX()

Notes:

• The Jacobian is defined as the PETSc type MatShell (https://petsc.org/release/docs/manualpages/Mat/MATSHELL.html)
• 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	Time stepping object (see PETSc TS)
t	Current time
Y	Solution vector
Ydot	Time-derivative of solution vector
a	Shift
A	Jacobian matrix
B	Preconditioning matrix
ctxt	Application context

Definition at line 61 of file PetscIJacobianIMEX.cpp.

{
  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);
}

PetscErrorCode PetscJacobianFunctionIMEX_JFNK	(	Mat	Jacobian,
		Vec	Y,
		Vec	F
	)

Computes the action of the Jacobian on a vector: See documentation for PetscIJacobianIMEX() for the definition of the Jacobian. This function computes its action on a vector for a nonlinear system by taking the directional derivative, as follows:

\begin{equation} {\bf f} = {\bf J}\left({\bf U}_0\right){\bf y} = \frac {\partial \mathcal{F}\left({\bf U}\right)} {\partial {\bf U}}\left({\bf U}_0\right) {\bf y} \approx \frac{1}{\epsilon} \left[ \mathcal{F}\left({\bf U}_0+\epsilon{\bf y}\right)-\mathcal{F}\left({\bf U}_0\right) \right] \end{equation}

In the context of the implicit-explicit (IMEX) time integration of the ODE given by

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

where \({\bf F}\) is the spatially discretized non-stiff right-hand-side (treated explicitly), \({\bf G}\) is the spatially discretized stiff right-hand-side (treated implicitly), we have that

\begin{equation} \mathcal{F}\left({\bf U}\right) \equiv \dot{\bf U} - {\bf G}\left({\bf U}\right) \Rightarrow {\bf J}\left({\bf U}\right) = \left[\alpha{\bf I} - \frac {\partial {\bf G}\left({\bf U}\right)} {\partial {\bf U}} \right]. \end{equation}

where \(\alpha\) (PETScContext::shift) is the shift variable specific to the time integration method. So this function computes

\begin{equation} {\bf f} = \alpha {\bf y} - \frac{1}{\epsilon} \left[ {\bf G}\left({\bf U}_0+\epsilon{\bf y}\right)-{\bf G}\left({\bf U}_0\right) \right] \end{equation}

In the code, \({\bf y}\) is Y, \(\bf f\) is F, and \({\bf U}_0\) is HyPar::uref (the reference solution at which the nonlinear Jacobian is computed). See papers on Jacobian-free Newton-Krylov (JFNK) methods to understand how \(\epsilon\) is computed.

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).

Notes:

• For a nonlinear spatial discretization, the right-hand-side must be computed without the nonlinearity (i.e. with a previously computed or "frozen" discretization operator). This ensures that the Jacobian being computed is consistent.
• 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

Jacobian	Jacobian matrix
Y	Input vector
F	Output vector (Jacobian times input vector

Definition at line 128 of file PetscIJacobianIMEX.cpp.

{
  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);
}

PetscErrorCode PetscJacobianFunctionIMEX_Linear	(	Mat	Jacobian,
		Vec	Y,
		Vec	F
	)

Computes the action of the Jacobian on a vector: See documentation for PetscIJacobianIMEX() for the definition of the Jacobian. This function computes its action on a vector for a linear system by taking the directional derivative, as follows:

\begin{equation} {\bf f} = {\bf J}{\bf y} = \frac {\partial \mathcal{F}\left({\bf U}\right)} {\partial {\bf U}} {\bf y} = \left[ \mathcal{F}\left({\bf U}_0+{\bf y}\right)-\mathcal{F}\left({\bf U}_0\right) \right], \end{equation}

where \(\mathcal{F}\) is linear, and thus \({\bf J}\) is a constant ( \(\mathcal{F}\left({\bf y}\right) = {\bf J}{\bf y}\)). In the context of the implicit-explicit (IMEX) time integration of the ODE given by

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

where \({\bf F}\) is the spatially discretized non-stiff right-hand-side (treated explicitly), \({\bf G}\) is the spatially discretized stiff right-hand-side (treated implicitly) and is linear, we have that

\begin{equation} \mathcal{F}\left({\bf U}\right) \equiv \dot{\bf U} - {\bf G}\left({\bf U}\right) \Rightarrow {\bf J}\left({\bf U}\right) = \left[\alpha{\bf I} - \frac {\partial {\bf G}\left({\bf U}\right)} {\partial {\bf U}} \right]. \end{equation}

where \(\alpha\) (PETScContext::shift) is the shift variable specific to the time integration method. So this function computes

\begin{equation} {\bf f} = \alpha {\bf y} - \left[ {\bf G}\left({\bf U}_0+{\bf y}\right)-{\bf G}\left({\bf U}_0\right) \right] \end{equation}

In the code, \({\bf y}\) is Y, \(\bf f\) is F, and \({\bf U}_0\) is HyPar::uref (the reference solution at which the nonlinear Jacobian is computed).

Since \(\mathcal{F}\) is linear,

\begin{equation} {\bf J}{\bf y} = \left[ \mathcal{F}\left({\bf U}_0+{\bf y}\right)-\mathcal{F}\left({\bf U}_0\right) \right] = \mathcal{F}\left({\bf y}\right). \end{equation}

However, the Jacobian is not computed as \(\mathcal{F}\left({\bf y}\right)\) because of the following reason: this function is used by the equation solver within the implicit time integrator in PETSc, and \({\bf y}\) represents the change in the solution, i.e. \(\Delta {\bf U}\), and not the solution \(\bf U\). Thus, evaluating \(\mathcal{F}\left({\bf y}\right)\) using HyPar::HyperbolicFunction, HyPar::ParabolicFunction, and HyPar::SourceFunction is incorrect since these functions expect \(\bf U\) as the input.

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).

Notes:

• For a nonlinear spatial discretization, the right-hand-side must be computed without the nonlinearity (i.e. with a previously computed or "frozen" discretization operator). This ensures that the Jacobian being computed is consistent, and is truly linear.
• 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

Jacobian	Jacobian matrix
Y	Input vector
F	Output vector (Jacobian times input vector

Definition at line 290 of file PetscIJacobianIMEX.cpp.

{
  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);
}