HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
petscinterface.h File Reference

Go to the source code of this file.

Functions

int TransferVecToPETSc (const double *const, Vec, void *, const int, const int)
 
int TransferVecFromPETSc (double *const, const Vec, void *, const int, const int)
 
int TransferMatToPETSc (void *, Mat, void *)
 
int PetscRegisterTIMethods (int)
 
PetscErrorCode PetscRHSFunctionExpl (TS, PetscReal, Vec, Vec, void *)
 
PetscErrorCode PetscRHSFunctionIMEX (TS, PetscReal, Vec, Vec, void *)
 
PetscErrorCode PetscIFunctionIMEX (TS, PetscReal, Vec, Vec, Vec, void *)
 
PetscErrorCode PetscIFunctionImpl (TS, PetscReal, Vec, Vec, Vec, void *)
 
PetscErrorCode PetscIJacobianIMEX (TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *)
 
PetscErrorCode PetscJacobianFunctionIMEX_JFNK (Mat, Vec, Vec)
 
PetscErrorCode PetscJacobianFunctionIMEX_Linear (Mat, Vec, Vec)
 
PetscErrorCode PetscIJacobian (TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *)
 
PetscErrorCode PetscJacobianFunction_JFNK (Mat, Vec, Vec)
 
PetscErrorCode PetscJacobianFunction_Linear (Mat, Vec, Vec)
 
int PetscGlobalDOF (void *)
 
int PetscCleanup (void *)
 
int PetscCreatePointList (void *)
 
PetscErrorCode PetscPreStage (TS, PetscReal)
 
PetscErrorCode PetscPostStage (TS, PetscReal, PetscInt, Vec *)
 
PetscErrorCode PetscPreTimeStep (TS)
 
PetscErrorCode PetscPostTimeStep (TS)
 
int PetscComputePreconMatIMEX (Mat, Vec, void *)
 
int PetscComputePreconMatImpl (Mat, Vec, void *)
 
int PetscJacobianMatNonzeroEntriesImpl (Mat, int, void *)
 
PetscErrorCode PetscTimeError (TS)
 
PetscErrorCode PetscSetInitialGuessROM (SNES, Vec, void *)
 

Function Documentation

int TransferVecToPETSc ( const double *const  u,
Vec  Y,
void *  ctxt,
const int  sim_idx,
const int  offset 
)

Copy data to a PETSc vector (used by PETSc time integrators, and with no ghost points) from a HyPar::u array (with ghost points).

See Also
TransferVecFromPETSc()
Parameters
uHyPar::u type array (with ghost points)
YPETSc vector
ctxtObject of type PETScContext
sim_idxSimulation object index
offsetOffset

Definition at line 24 of file TransferToPETSc.cpp.

29 {
30  PETScContext* context = (PETScContext*) ctxt;
31  SimulationObject* sim = (SimulationObject*) context->simobj;
32  double* Yarr;
33 
34  PetscFunctionBegin;
35  VecGetArray(Y,&Yarr);
36  std::vector<int> index(sim[sim_idx].solver.ndims,0);
37  ArrayCopynD( sim[sim_idx].solver.ndims,
38  u,
39  (Yarr+offset),
40  sim[sim_idx].solver.dim_local,
41  sim[sim_idx].solver.ghosts,
42  0,
43  index.data(),
44  sim[sim_idx].solver.nvars );
45  VecRestoreArray(Y,&Yarr);
46 
47  PetscFunctionReturn(0);
48 }
INLINE int ArrayCopynD(int, const double *, double *, int *, int, int, int *, int)
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int TransferVecFromPETSc ( double *const  u,
const Vec  Y,
void *  ctxt,
const int  sim_idx,
const int  offset 
)

Copy data from a PETSc vector (used by PETSc time integrators, and with no ghost points) to a HyPar::u array (with ghost points).

See Also
TransferVecToPETSc()
Parameters
uHyPar::u type array (with ghost points)
YPETSc vector
ctxtObject of type PETScContext
sim_idxSimulation object index
offsetOffset

Definition at line 23 of file TransferFromPETSc.cpp.

28 {
29  PETScContext* context = (PETScContext*) ctxt;
30  SimulationObject* sim = (SimulationObject*) context->simobj;
31  const double* Yarr;
32 
33  PetscFunctionBegin;
34  VecGetArrayRead(Y,&Yarr);
35  std::vector<int> index(sim[sim_idx].solver.ndims,0);
36  ArrayCopynD( sim[sim_idx].solver.ndims,
37  (Yarr+offset),
38  u,
39  sim[sim_idx].solver.dim_local,
40  0,
41  sim[sim_idx].solver.ghosts,
42  index.data(),
43  sim[sim_idx].solver.nvars );
44  VecRestoreArrayRead(Y,&Yarr);
45 
46  PetscFunctionReturn(0);
47 }
INLINE int ArrayCopynD(int, const double *, double *, int *, int, int, int *, int)
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int TransferMatToPETSc ( void *  J,
Mat  A,
void *  ctxt 
)

Copy a matrix of type BandedMatrix to a PETSc matrix.

Parameters
JMatrix of type BandedMatrix
APETSc matrix
ctxtObject of type PETScContext

Definition at line 53 of file TransferToPETSc.cpp.

56 {
57  BandedMatrix *M = (BandedMatrix*) J;
58  PetscErrorCode ierr = 0;
59  int bs = M->BlockSize, nbands = M->nbands, bs2 = bs*bs;
60 
61  for (int i=0; i<M->nrows_local; i++) {
62  int colind[nbands];
63  double val[bs][bs*nbands];
64  for (int n=0; n<nbands; n++) {
65  colind[n] = M->ncol[nbands*i+n];
66  for (int p=0; p<bs; p++) {
67  for (int q = 0; q<bs; q++) {
68  val[p][n*bs+q] = M->data[i*nbands*bs2+n*bs2+p*bs+q];
69  }
70  }
71  }
72  MatSetValuesBlocked(A,1,&M->nrow[i],M->nbands,&colind[0],&val[0][0],INSERT_VALUES);
73  }
74 
75  MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
76  MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
77 
78  return(0);
79 }
double * data
Definition: bandedmatrix.h:24
Structure for defining a banded block matrix.
Definition: bandedmatrix.h:18
int PetscRegisterTIMethods ( int  rank)

This function allows the registering of custom multi-stage time integration methods and using it with PETSc. More than one method may be provided, each with a unique name (the name must not conflict with the names of existing methods in PETSc). The methods should be provided in a file "time_method.inp", and the following must be provided:

  • Name
  • Class (i.e. RK or ARK)
  • Number of stages
  • Theoretical order
  • Butcher tableaux

The method can then be used by invoking it by its name. For example, if a custom ARKIMEX method "foo" is provided through "time_method.inp", it can be used by specifying "-ts_type arkimex -ts_arkimex_type foo" in the command line or ".petscrc" file.

See /Examples/PETScInputs for examples of the input files required to provide a time integration method.

Currently supported classes of time integrators for which custom methods can be provided:

To do:

  • Add support for TSGLEE when it gets merged to PETSc's master.
Parameters
rankMPI rank

Definition at line 46 of file PetscRegisterTIMethods.cpp.

47 {
48  PetscErrorCode ierr;
49  int ierr2;
50 
51  PetscFunctionBegin;
52 
53  /* Note: all processors read and register the custom methods */
54  /* instead of root doing it and broadcasting. */
55  FILE *in;
56  in = fopen("time_method.inp","r");
57  if (in) {
58  int n, N; /* N = total number of methods specified in the file */
59  ierr2 = fscanf(in,"%d",&N); if (ierr2 != 1) return(1);
60  for (n = 0; n < N; n++) {
61  char name[_MAX_STRING_SIZE_]; /* name of the scheme */
62  char type[_MAX_STRING_SIZE_]; /* type of scheme - ARKIMEX or RK */
63  PetscInt s, order; /* number of stages and order */
64  PetscReal *A , *b, *c; /* Butcher tableaux entries for non-stiff terms */
65  PetscReal *At, *bt, *ct; /* Butcher tableaux entries for the stiff terms */
66  PetscReal *bemb, *bembt; /* Embedded method coefficients */
67  PetscInt pinterp; /* order of interpolation scheme */
68  PetscReal *bint,*bintt; /* Dense output interpolation coefficients */
69 
70  /* Initializations */
71  strcpy(name,"");
72  s = order = pinterp = 0;
73  A = b = c = NULL;
74  At = bt = ct = NULL;
75  bemb = bembt = NULL;
76  bint = bintt = NULL;
77 
78  /* Read the method */
79  char word[_MAX_STRING_SIZE_];
80  ierr2 = fscanf(in,"%s",word); if (ierr2 != 1) return(1);
81  if (!strcmp(word,"begin")) {
82  while (strcmp(word, "end")) {
83  ierr2 = fscanf(in,"%s",word); if (ierr2 != 1) return(1);
84  if (!strcmp(word,"name")) { ierr2 = fscanf(in,"%s",name); if (ierr2 != 1) return(1); }
85  else if (!strcmp(word,"class")) { ierr2 = fscanf(in,"%s",type); if (ierr2 != 1) return(1); }
86  else if (!strcmp(word,"nstages")) { ierr2 = fscanf(in,"%d",&s); if (ierr2 != 1) return(1); }
87  else if (!strcmp(word,"order")) { ierr2 = fscanf(in,"%d",&order); if (ierr2 != 1) return(1); }
88  else if (!strcmp(word,"pinterp")) { ierr2 = fscanf(in,"%d",&pinterp);if (ierr2 != 1) return(1); }
89  else if (!strcmp(word,"At")) {
90  if (s == 0) {
91  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
92  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
93  return(1);
94  } else {
95  At = (PetscReal*) calloc (s*s, sizeof(PetscReal));
96  int i, j;
97  for (i = 0; i < s; i++) {
98  for (j = 0; j < s; j++) {
99  ierr2 = fscanf(in,"%lf",&At[i*s+j]); if (ierr2 != 1) return(1);
100  }
101  }
102  }
103  } else if (!strcmp(word,"A")) {
104  if (s == 0) {
105  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
106  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
107  return(1);
108  } else {
109  A = (PetscReal*) calloc (s*s, sizeof(PetscReal));
110  int i, j;
111  for (i = 0; i < s; i++) {
112  for (j = 0; j < s; j++) {
113  ierr2 = fscanf(in,"%lf",&A[i*s+j]); if (ierr2 != 1) return(1);
114  }
115  }
116  }
117  } else if (!strcmp(word,"bt")) {
118  if (s == 0) {
119  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
120  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
121  return(1);
122  } else {
123  bt = (PetscReal*) calloc (s, sizeof(PetscReal));
124  int i;
125  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&bt[i]); if (ierr2 != 1) return(1);
126  }
127  } else if (!strcmp(word,"b")) {
128  if (s == 0) {
129  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
130  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
131  return(1);
132  } else {
133  b = (PetscReal*) calloc (s, sizeof(PetscReal));
134  int i;
135  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&b[i]); if (ierr2 != 1) return(1);
136  }
137  } else if (!strcmp(word,"ct")) {
138  if (s == 0) {
139  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
140  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
141  return(1);
142  } else {
143  ct = (PetscReal*) calloc (s, sizeof(PetscReal));
144  int i;
145  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&ct[i]); if (ierr2 != 1) return(1);
146  }
147  } else if (!strcmp(word,"c")) {
148  if (s == 0) {
149  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
150  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
151  return(1);
152  } else {
153  c = (PetscReal*) calloc (s, sizeof(PetscReal));
154  int i;
155  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&c[i]); if (ierr2 != 1) return(1);
156  }
157  } else if (!strcmp(word,"bembt")) {
158  if (s == 0) {
159  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
160  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
161  return(1);
162  } else {
163  bembt = (PetscReal*) calloc (s, sizeof(PetscReal));
164  int i;
165  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&bembt[i]); if (ierr2 != 1) return(1);
166  }
167  } else if (!strcmp(word,"bemb")) {
168  if (s == 0) {
169  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages must be defined ");
170  if (!rank) fprintf(stderr,"before specifying the Butcher tableaux entries.\n" );
171  return(1);
172  } else {
173  bemb = (PetscReal*) calloc (s, sizeof(PetscReal));
174  int i;
175  for (i = 0; i < s; i++) ierr2 = fscanf(in,"%lf",&bemb[i]); if (ierr2 != 1) return(1);
176  }
177  } else if (!strcmp(word,"bintt")) {
178  if (s == 0 || pinterp == 0) {
179  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages and pinterp must be " );
180  if (!rank) fprintf(stderr,"defined as positive values before specifying interpolation coeffs.\n");
181  return(1);
182  } else {
183  bintt = (PetscReal*) calloc (s*pinterp, sizeof(PetscReal));
184  int i, j;
185  for (i = 0; i < s; i++) {
186  for (j = 0; j < pinterp; j++) {
187  ierr2 = fscanf(in,"%lf",&bintt[i*s+j]); if (ierr2 != 1) return(1);
188  }
189  }
190  }
191  } else if (!strcmp(word,"bint")) {
192  if (s == 0 || pinterp == 0) {
193  if (!rank) fprintf(stderr,"Error in PetscRegisterTIMethods(): nstages and pinterp must be " );
194  if (!rank) fprintf(stderr,"defined as positive values before specifying interpolation coeffs.\n");
195  return(1);
196  } else {
197  bint = (PetscReal*) calloc (s*pinterp, sizeof(PetscReal));
198  int i, j;
199  for (i = 0; i < s; i++) {
200  for (j = 0; j < pinterp; j++) {
201  ierr2 = fscanf(in,"%lf",&bint[i*s+j]); if (ierr2 != 1) return(1);
202  }
203  }
204  }
205  }
206  }
207  } else {
208  if (!rank) fprintf(stderr,"Error: Illegal format in file \"time_method.inp\" (expected keyword \"begin\").\n");
209  return(1);
210  }
211 
212  /* Register the method */
213  if (!strcmp(type,"arkimex")) {
214  if (A && At) {
215  ierr = TSARKIMEXRegister(name,order,s,At,bt,ct,A,b,c,bembt,bemb,pinterp,bintt,bint); CHKERRQ(ierr);
216  if (!rank) {
217  printf("\nRegistered custom ARKIMEX scheme \"%s\" with the following Butcher tableaux:-\n",name);
218  int i,j;
219  for (i = 0; i < s; i++) {
220  if (c) printf(" %+1.5lf |",c[i]);
221  else printf(" |");
222  for (j = 0; j < s; j++) printf (" %+1.5lf :",A[i*s+j]);
223  printf(" ");
224  if (ct) printf("%+1.5lf |",ct[i]);
225  else printf(" |");
226  for (j = 0; j < s; j++) printf (" %+1.5lf :",At[i*s+j]);
227  printf("\n");
228  }
229  printf(" ---------|");
230  for (j = 0; j < s; j++) printf("-----------");
231  printf(" ");
232  printf("---------|");
233  for (j = 0; j < s; j++) printf("-----------");
234  printf("\n");
235  printf(" |");
236  if (b) for (j = 0; j < s; j++) printf(" %+1.5lf :",b[j]);
237  else for (j = 0; j < s; j++) printf(" :");
238  printf(" ");
239  printf(" |");
240  if (bt) for (j = 0; j < s; j++) printf(" %+1.5lf :",bt[j]);
241  else for (j = 0; j < s; j++) printf(" :");
242  printf("\n\n");
243  }
244  } else {
245  if (!rank) {
246  fprintf(stderr,"Warning in PetscRegisterTIMethods(): Failed to register method ");
247  fprintf(stderr,"(A or At not defined).\n");
248  }
249  }
250  } else if (!strcmp(type,"rk")) {
251  if (A) {
252  ierr = TSRKRegister(name,order,s,A,b,c,bemb,pinterp,bint); CHKERRQ(ierr);
253  if (!rank) {
254  printf("\nRegistered custom RK scheme \"%s\" with the following Butcher tableaux:-\n",name);
255  int i,j;
256  for (i = 0; i < s; i++) {
257  if (c) printf(" %+1.5lf |",c[i]);
258  else printf(" |");
259  for (j = 0; j < s; j++) printf (" %+1.5lf :",A[i*s+j]);
260  printf("\n");
261  }
262  printf(" ---------|");
263  for (j = 0; j < s; j++) printf("-----------");
264  printf("\n");
265  printf(" |");
266  if (b) for (j = 0; j < s; j++) printf(" %+1.5lf :",b[j]);
267  else for (j = 0; j < s; j++) printf(" :");
268  printf("\n\n");
269  }
270  } else {
271  if (!rank) {
272  fprintf(stderr,"Warning in PetscRegisterTIMethods(): Failed to register method ");
273  fprintf(stderr,"(A not defined).\n");
274  }
275  }
276  } else {
277  if (!rank){
278  fprintf(stderr,"Error in PetscRegisterTIMethods(): %s class of time-integration schemes ",type);
279  fprintf(stderr,"does not support custom method registration and usage.\n");
280  }
281  }
282 
283  /* Free the arrays */
284  if (At) free(At);
285  if (bt) free(bt);
286  if (ct) free(ct);
287  if (A) free(A);
288  if (b) free(b);
289  if (c) free(c);
290  if (bembt) free(bembt);
291  if (bemb) free(bemb);
292  if (bintt) free(bintt);
293  if (bint) free(bint);
294 
295  }
296  }
297  PetscFunctionReturn(0);
298 }
#define _MAX_STRING_SIZE_
Definition: basic.h:14
PetscErrorCode PetscRHSFunctionExpl ( TS  ts,
PetscReal  t,
Vec  Y,
Vec  F,
void *  ctxt 
)

Compute the right-hand-side (RHS) for the explicit time integration of the governing equations: The spatially discretized ODE can be expressed as

\begin{equation} \frac {d{\bf U}} {dt} = {\bf F}\left({\bf U}\right). \end{equation}

This function computes \({\bf F}\left({\bf U}\right)\), given \({\bf U}\).

See Also
TSSetRHSFunction (https://petsc.org/release/docs/manualpages/TS/TSSetRHSFunction.html)

Note:

Parameters
tsTime integration object
tCurrent simulation time
YState vector (input)
FThe computed right-hand-side vector
ctxtObject of type PETScContext

Definition at line 36 of file PetscRHSFunctionExpl.cpp.

41 {
42  PETScContext* context = (PETScContext*) ctxt;
43  SimulationObject* sim = (SimulationObject*) context->simobj;
44  int nsims = context->nsims;
45 
46  PetscFunctionBegin;
47 
48  for (int ns = 0; ns < nsims; ns++) {
49 
50  HyPar* solver = &(sim[ns].solver);
51  MPIVariables* mpi = &(sim[ns].mpi);
52 
53  solver->count_RHSFunction++;
54 
55  int size = solver->npoints_local_wghosts;
56  double* u = solver->u;
57  double* rhs = solver->rhs;
58 
59  /* copy solution from PETSc vector */
60  TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
61  /* apply boundary conditions and exchange data over MPI interfaces */
62  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
64  solver->nvars,
65  solver->dim_local,
66  solver->ghosts,
67  mpi,
68  u );
69 
70  /* Evaluate hyperbolic, parabolic and source terms and the RHS */
71  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,1,solver->FFunction,solver->Upwind);
72 
73  solver->ParabolicFunction (solver->par,u,solver,mpi,t);
74  solver->SourceFunction (solver->source,u,solver,mpi,t);
75 
76  _ArraySetValue_(rhs,size*solver->nvars,0.0);
77  _ArrayAXPY_(solver->hyp ,-1.0,rhs,size*solver->nvars);
78  _ArrayAXPY_(solver->par , 1.0,rhs,size*solver->nvars);
79  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
80 
81  /* Transfer RHS to PETSc vector */
82  TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
83 
84  }
85 
86  PetscFunctionReturn(0);
87 }
double * hyp
Definition: hypar.h:119
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int * dim_local
Definition: hypar.h:37
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
int count_RHSFunction
Definition: hypar.h:422
double * par
Definition: hypar.h:122
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399
PetscErrorCode PetscRHSFunctionIMEX ( TS  ts,
PetscReal  t,
Vec  Y,
Vec  F,
void *  ctxt 
)

Compute the explicitly-treated right-hand-side (RHS) 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 F}\left({\bf U}\right)\) represents all terms that the user has indicated to be integrated in time explicitly (PETScContext::flag_hyperbolic_f, PETScContext::flag_hyperbolic_df, PETScContext::flag_hyperbolic, PETScContext::flag_parabolic, and PETScContext::flag_source).

This function computes \({\bf F}\left({\bf U}\right)\), given \({\bf U}\).

See Also
PetscIFunctionIMEX(), TSSetRHSFunction (https://petsc.org/release/docs/manualpages/TS/TSSetRHSFunction.html)

Note:

Parameters
tsTime integration object
tCurrent simulation time
YState vector (input)
FThe computed right-hand-side vector
ctxtObject of type PETScContext

Definition at line 49 of file PetscRHSFunctionIMEX.cpp.

54 {
55  PETScContext* context = (PETScContext*) ctxt;
56  SimulationObject* sim = (SimulationObject*) context->simobj;
57  int nsims = context->nsims;
58 
59  context->waqt = t;
60 
61  PetscFunctionBegin;
62  for (int ns = 0; ns < nsims; ns++) {
63 
64  HyPar* solver = &(sim[ns].solver);
65  MPIVariables* mpi = &(sim[ns].mpi);
66 
67  solver->count_RHSFunction++;
68 
69  int size = solver->npoints_local_wghosts;
70  double *u = solver->u;
71  double *rhs = solver->rhs;
72 
73  /* copy solution from PETSc vector */
74  TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
75  /* apply boundary conditions and exchange data over MPI interfaces */
76  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
77  MPIExchangeBoundariesnD(solver->ndims,solver->nvars,solver->dim_local,
78  solver->ghosts,mpi,u);
79 
80  /* initialize right-hand side to zero */
81  _ArraySetValue_(rhs,size*solver->nvars,0.0);
82 
83  /* Evaluate hyperbolic, parabolic and source terms and the RHS */
84  if ((!strcmp(solver->SplitHyperbolicFlux,"yes")) && solver->flag_fdf_specified) {
85  if (context->flag_hyperbolic_f == _EXPLICIT_) {
86  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FdFFunction,solver->UpwindFdF);
87  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
88  }
89  if (context->flag_hyperbolic_df == _EXPLICIT_) {
90  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF);
91  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
92  }
93  } else if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
94  if (context->flag_hyperbolic_f == _EXPLICIT_) {
95  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FFunction,solver->Upwind);
96  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
97  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF);
98  _ArrayAXPY_(solver->hyp, 1.0,rhs,size*solver->nvars);
99  }
100  if (context->flag_hyperbolic_df == _EXPLICIT_) {
101  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF);
102  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
103  }
104  } else {
105  if (context->flag_hyperbolic == _EXPLICIT_) {
106  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FFunction,solver->Upwind);
107  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
108  }
109  }
110  if (context->flag_parabolic == _EXPLICIT_) {
111  solver->ParabolicFunction (solver->par,u,solver,mpi,t);
112  _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
113  }
114  if (context->flag_source == _EXPLICIT_) {
115  solver->SourceFunction (solver->source,u,solver,mpi,t);
116  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
117  }
118 
119  /* Transfer RHS to PETSc vector */
120  TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
121  }
122 
123  PetscFunctionReturn(0);
124 }
double * hyp
Definition: hypar.h:119
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
int(* dFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:280
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
#define _EXPLICIT_
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int flag_fdf_specified
Definition: hypar.h:290
int(* UpwinddF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:300
int * dim_local
Definition: hypar.h:37
int(* UpwindFdF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:307
int(* FdFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:286
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
int count_RHSFunction
Definition: hypar.h:422
double * par
Definition: hypar.h:122
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399
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:

Parameters
tsThe time integration object
tCurrent solution time
YState vector (input)
YdotTime derivative of the state vector (input)
FThe computed function vector
ctxtObject of type PETScContext

Definition at line 53 of file PetscIFunctionIMEX.cpp.

59 {
60  PETScContext* context = (PETScContext*) ctxt;
61  SimulationObject* sim = (SimulationObject*) context->simobj;
62  int nsims = context->nsims;
63 
64  PetscFunctionBegin;
65 
66  context->waqt = t;
67 
68  for (int ns = 0; ns < nsims; ns++) {
69 
70  HyPar* solver = &(sim[ns].solver);
71  MPIVariables* mpi = &(sim[ns].mpi);
72 
73  solver->count_IFunction++;
74 
75  int size = solver->npoints_local_wghosts;
76  double *u = solver->u;
77  double *rhs = solver->rhs;
78 
79  /* copy solution from PETSc vector */
80  TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
81  /* apply boundary conditions and exchange data over MPI interfaces */
82  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
84  solver->nvars,
85  solver->dim_local,
86  solver->ghosts,
87  mpi,
88  u );
89 
90  /* initialize right-hand side to zero */
91  _ArraySetValue_(rhs,size*solver->nvars,0.0);
92 
93  /* Evaluate hyperbolic, parabolic and source terms and the RHS */
94  if ((!strcmp(solver->SplitHyperbolicFlux,"yes")) && solver->flag_fdf_specified) {
95  if (context->flag_hyperbolic_f == _IMPLICIT_) {
96  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FdFFunction,solver->UpwindFdF);
97  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
98  }
99  if (context->flag_hyperbolic_df == _IMPLICIT_) {
100  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF);
101  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
102  }
103  } else if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
104  if (context->flag_hyperbolic_f == _IMPLICIT_) {
105  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FFunction,solver->Upwind);
106  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
107  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF);
108  _ArrayAXPY_(solver->hyp, 1.0,rhs,size*solver->nvars);
109  }
110  if (context->flag_hyperbolic_df == _IMPLICIT_) {
111  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->dFFunction,solver->UpwinddF);
112  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
113  }
114  } else {
115  if (context->flag_hyperbolic == _IMPLICIT_) {
116  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,0,solver->FFunction,solver->Upwind);
117  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
118  }
119  }
120  if (context->flag_parabolic == _IMPLICIT_) {
121  solver->ParabolicFunction (solver->par,u,solver,mpi,t);
122  _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
123  }
124  if (context->flag_source == _IMPLICIT_) {
125  solver->SourceFunction (solver->source,u,solver,mpi,t);
126  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
127  }
128 
129  /* save a copy of the solution and RHS for use in IJacobian */
130  _ArrayCopy1D_(u ,solver->uref ,(size*solver->nvars));
131  _ArrayCopy1D_(rhs,solver->rhsref,(size*solver->nvars));
132 
133  /* Transfer RHS to PETSc vector */
134  TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
135 
136  }
137 
138  /* LHS = Ydot - F(u) */
139  VecAYPX(F,-1.0,Ydot);
140 
141  PetscFunctionReturn(0);
142 }
double * hyp
Definition: hypar.h:119
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
int(* dFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:280
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int flag_fdf_specified
Definition: hypar.h:290
int(* UpwinddF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:300
int * dim_local
Definition: hypar.h:37
int(* UpwindFdF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:307
int(* FdFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:286
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
#define _IMPLICIT_
double * par
Definition: hypar.h:122
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
double * rhsref
Definition: hypar.h:398
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
double * uref
Definition: hypar.h:397
#define _ArrayCopy1D_(x, y, size)
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int nvars
Definition: hypar.h:29
int count_IFunction
Definition: hypar.h:422
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399
PetscErrorCode PetscIFunctionImpl ( TS  ts,
PetscReal  t,
Vec  Y,
Vec  Ydot,
Vec  F,
void *  ctxt 
)

Compute the left-hand-side for the implicit time integration of the governing equations: The spatially discretized ODE can be expressed as

\begin{equation} \frac {d{\bf U}} {dt} = {\bf F}\left({\bf U}\right). \end{equation}

This function computes \(\dot{\bf U} - {\bf F}\left({\bf U}\right)\), given \({\bf U},\dot{\bf U}\).

See Also
TSSetIFunction (https://petsc.org/release/docs/manualpages/TS/TSSetIFunction.html)

Note:

Parameters
tsTime integration object
tCurrent simulation time
YState vector (input)
YdotTime derivative of the state vector (input)
FThe computed right-hand-side vector
ctxtObject of type PETScContext

Definition at line 37 of file PetscIFunctionImpl.cpp.

43 {
44  PETScContext* context = (PETScContext*) ctxt;
45  SimulationObject* sim = (SimulationObject*) context->simobj;
46  int nsims = context->nsims;
47 
48  PetscFunctionBegin;
49  for (int ns = 0; ns < nsims; ns++) {
50 
51  HyPar* solver = &(sim[ns].solver);
52  MPIVariables* mpi = &(sim[ns].mpi);
53 
54  solver->count_RHSFunction++;
55 
56  int size = solver->npoints_local_wghosts;
57  double* u = solver->u;
58  double* rhs = solver->rhs;
59 
60  /* copy solution from PETSc vector */
61  TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
62  /* apply boundary conditions and exchange data over MPI interfaces */
63  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
65  solver->nvars,
66  solver->dim_local,
67  solver->ghosts,
68  mpi,
69  u );
70 
71  /* Evaluate hyperbolic, parabolic and source terms and the RHS */
72  solver->HyperbolicFunction(solver->hyp,u,solver,mpi,t,1,solver->FFunction,solver->Upwind);
73  solver->ParabolicFunction (solver->par,u,solver,mpi,t);
74  solver->SourceFunction (solver->source,u,solver,mpi,t);
75 
76  _ArraySetValue_(rhs,size*solver->nvars,0.0);
77  _ArrayAXPY_(solver->hyp ,-1.0,rhs,size*solver->nvars);
78  _ArrayAXPY_(solver->par , 1.0,rhs,size*solver->nvars);
79  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
80 
81  /* save a copy of the solution and RHS for use in IJacobian */
82  _ArrayCopy1D_(u ,solver->uref ,(size*solver->nvars));
83  _ArrayCopy1D_(rhs,solver->rhsref,(size*solver->nvars));
84 
85  /* Transfer RHS to PETSc vector */
86  TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
87 
88  }
89 
90  /* LHS = Ydot - F(u) */
91  VecAYPX(F,-1.0,Ydot);
92 
93  PetscFunctionReturn(0);
94 }
double * hyp
Definition: hypar.h:119
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int * dim_local
Definition: hypar.h:37
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
int count_RHSFunction
Definition: hypar.h:422
double * par
Definition: hypar.h:122
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
double * rhsref
Definition: hypar.h:398
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
double * uref
Definition: hypar.h:397
#define _ArrayCopy1D_(x, y, size)
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399
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:

See Also
PetscIFunctionIMEX()

Notes:

Parameters
tsTime stepping object (see PETSc TS)
tCurrent time
YSolution vector
YdotTime-derivative of solution vector
aShift
AJacobian matrix
BPreconditioning matrix
ctxtApplication context

Definition at line 61 of file PetscIJacobianIMEX.cpp.

69 {
70  PETScContext* context = (PETScContext*) ctxt;
71 
72  PetscFunctionBegin;
73  for (int ns = 0; ns < context->nsims; ns++) {
74  ((SimulationObject*)context->simobj)[ns].solver.count_IJacobian++;
75  }
76  context->shift = a;
77  context->waqt = t;
78  /* Construct preconditioning matrix */
79  if (context->flag_use_precon) PetscComputePreconMatIMEX(B,Y,context);
80 
81  PetscFunctionReturn(0);
82 }
int PetscComputePreconMatIMEX(Mat, Vec, void *)
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
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
JacobianJacobian matrix
YInput vector
FOutput vector (Jacobian times input vector

Definition at line 128 of file PetscIJacobianIMEX.cpp.

131 {
132  PETScContext* context(nullptr);
133 
134  PetscFunctionBegin;
135 
136  MatShellGetContext(Jacobian,&context);
137  SimulationObject* sim = (SimulationObject*) context->simobj;
138  int nsims = context->nsims;
139 
140  double normY;
141  VecNorm(Y,NORM_2,&normY);
142 
143  if (normY < 1e-16) {
144 
145  /* F = 0 */
146  VecZeroEntries(F);
147  /* [J]Y = aY - F(Y) */
148  VecAXPBY(F,context->shift,0,Y);
149 
150  } else {
151 
152  double epsilon = context->jfnk_eps / normY;
153  double t = context->waqt; /* current stage/step time */
154 
155  for (int ns = 0; ns < nsims; ns++) {
156 
157  HyPar* solver = &(sim[ns].solver);
158  MPIVariables* mpi = &(sim[ns].mpi);
159  solver->count_IJacFunction++;
160 
161  int size = solver->npoints_local_wghosts;
162 
163  double *u = solver->u;
164  double *uref = solver->uref;
165  double *rhsref = solver->rhsref;
166  double *rhs = solver->rhs;
167 
168  /* copy solution from PETSc vector */
169  TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
170  _ArrayAYPX_(uref,epsilon,u,size*solver->nvars);
171  /* apply boundary conditions and exchange data over MPI interfaces */
172  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
174  solver->nvars,
175  solver->dim_local,
176  solver->ghosts,
177  mpi,
178  u );
179 
180  /* Evaluate hyperbolic, parabolic and source terms and the RHS for U+dU */
181  _ArraySetValue_(rhs,size*solver->nvars,0.0);
182  if ((!strcmp(solver->SplitHyperbolicFlux,"yes")) && solver->flag_fdf_specified) {
183  if (context->flag_hyperbolic_f == _IMPLICIT_) {
184  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
185  solver->FdFFunction,solver->UpwindFdF);
186  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
187  }
188  if (context->flag_hyperbolic_df == _IMPLICIT_) {
189  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
190  solver->dFFunction,solver->UpwinddF);
191  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
192  }
193  } else if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
194  if (context->flag_hyperbolic_f == _IMPLICIT_) {
195  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
196  solver->FFunction,solver->Upwind);
197  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
198  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
199  solver->dFFunction,solver->UpwinddF);
200  _ArrayAXPY_(solver->hyp, 1.0,rhs,size*solver->nvars);
201  }
202  if (context->flag_hyperbolic_df == _IMPLICIT_) {
203  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
204  solver->dFFunction,solver->UpwinddF);
205  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
206  }
207  } else {
208  if (context->flag_hyperbolic == _IMPLICIT_) {
209  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
210  solver->FFunction,solver->Upwind);
211  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
212  }
213  }
214  if (context->flag_parabolic == _IMPLICIT_) {
215  solver->ParabolicFunction (solver->par,u,solver,mpi,t);
216  _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
217  }
218  if (context->flag_source == _IMPLICIT_) {
219  solver->SourceFunction (solver->source,u,solver,mpi,t);
220  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
221  }
222 
223  _ArrayAXPY_(rhsref,-1.0,rhs,size*solver->nvars);
224  /* Transfer RHS to PETSc vector */
225  TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
226  }
227 
228  /* [J]Y = aY - F(Y) */
229  VecAXPBY(F,context->shift,(-1.0/epsilon),Y);
230 
231  }
232 
233  PetscFunctionReturn(0);
234 }
double * hyp
Definition: hypar.h:119
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
int(* dFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:280
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int flag_fdf_specified
Definition: hypar.h:290
int(* UpwinddF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:300
int * dim_local
Definition: hypar.h:37
int(* UpwindFdF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:307
int(* FdFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:286
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
#define _IMPLICIT_
double * par
Definition: hypar.h:122
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
double * rhsref
Definition: hypar.h:398
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
double * uref
Definition: hypar.h:397
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int count_IJacFunction
Definition: hypar.h:422
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
#define _ArrayAYPX_(x, a, y, size)
Structure defining a simulation.
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399
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
JacobianJacobian matrix
YInput vector
FOutput vector (Jacobian times input vector

Definition at line 290 of file PetscIJacobianIMEX.cpp.

294 {
295  PETScContext* context(nullptr);
296 
297  PetscFunctionBegin;
298 
299  MatShellGetContext(Jacobian,&context);
300  SimulationObject* sim = (SimulationObject*) context->simobj;
301  int nsims = context->nsims;
302 
303  double normY;
304  VecNorm(Y,NORM_2,&normY);
305 
306  if (normY < 1e-16) {
307 
308  /* F = 0 */
309  VecZeroEntries(F);
310  /* [J]Y = aY - F(Y) */
311  VecAXPBY(F,context->shift,0,Y);
312 
313  } else {
314 
315  double t = context->waqt; /* current stage/step time */
316 
317  for (int ns = 0; ns < nsims; ns++) {
318 
319  HyPar* solver = &(sim[ns].solver);
320  MPIVariables* mpi = &(sim[ns].mpi);
321  solver->count_IJacFunction++;
322 
323  int size = solver->npoints_local_wghosts;
324 
325  double *u = solver->u;
326  double *uref = solver->uref;
327  double *rhsref = solver->rhsref;
328  double *rhs = solver->rhs;
329 
330  /* copy solution from PETSc vector */
331  TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
332  _ArrayAYPX_(uref,1.0,u,size*solver->nvars);
333  /* apply boundary conditions and exchange data over MPI interfaces */
334  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
336  solver->nvars,
337  solver->dim_local,
338  solver->ghosts,
339  mpi,
340  u );
341 
342  /* Evaluate hyperbolic, parabolic and source terms and the RHS for U+dU */
343  _ArraySetValue_(rhs,size*solver->nvars,0.0);
344  if ((!strcmp(solver->SplitHyperbolicFlux,"yes")) && solver->flag_fdf_specified) {
345  if (context->flag_hyperbolic_f == _IMPLICIT_) {
346  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
347  solver->FdFFunction,solver->UpwindFdF);
348  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
349  }
350  if (context->flag_hyperbolic_df == _IMPLICIT_) {
351  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
352  solver->dFFunction,solver->UpwinddF);
353  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
354  }
355  } else if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
356  if (context->flag_hyperbolic_f == _IMPLICIT_) {
357  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
358  solver->FFunction,solver->Upwind);
359  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
360  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
361  solver->dFFunction,solver->UpwinddF);
362  _ArrayAXPY_(solver->hyp, 1.0,rhs,size*solver->nvars);
363  }
364  if (context->flag_hyperbolic_df == _IMPLICIT_) {
365  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
366  solver->dFFunction,solver->UpwinddF);
367  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
368  }
369  } else {
370  if (context->flag_hyperbolic == _IMPLICIT_) {
371  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
372  solver->FFunction,solver->Upwind);
373  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
374  }
375  }
376  if (context->flag_parabolic == _IMPLICIT_) {
377  solver->ParabolicFunction (solver->par,u,solver,mpi,t);
378  _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
379  }
380  if (context->flag_source == _IMPLICIT_) {
381  solver->SourceFunction (solver->source,u,solver,mpi,t);
382  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
383  }
384 
385  _ArrayAXPY_(rhsref,-1.0,rhs,size*solver->nvars);
386  /* Transfer RHS to PETSc vector */
387  TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
388  }
389 
390  /* [J]Y = aY - F(Y) */
391  VecAXPBY(F,context->shift,-1.0,Y);
392 
393  }
394 
395  PetscFunctionReturn(0);
396 }
double * hyp
Definition: hypar.h:119
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
int(* dFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:280
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int flag_fdf_specified
Definition: hypar.h:290
int(* UpwinddF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:300
int * dim_local
Definition: hypar.h:37
int(* UpwindFdF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:307
int(* FdFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:286
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
#define _IMPLICIT_
double * par
Definition: hypar.h:122
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
double * rhsref
Definition: hypar.h:398
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
double * uref
Definition: hypar.h:397
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int count_IJacFunction
Definition: hypar.h:422
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
#define _ArrayAYPX_(x, a, y, size)
Structure defining a simulation.
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399
PetscErrorCode PetscIJacobian ( TS  ts,
PetscReal  t,
Vec  Y,
Vec  Ydot,
PetscReal  a,
Mat  A,
Mat  B,
void *  ctxt 
)

Compute the Jacobian for implicit 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) \Rightarrow \dot{\bf U} - {\bf F}\left({\bf U}\right) = 0, \end{equation}

where \({\bf F}\) is the spatially discretized right-hand-side, and \({\bf U}\) represents the entire solution vector. The Jacobian is thus given by:

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

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

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 (PetscJacobianFunction_JFNK() for nonlinear systems, and PetscJacobianFunction_Linear() for linear systems). This approach works well with iterative solvers. Thus, this function just does the following:

Notes:

Parameters
tsTime stepping object (see PETSc TS)
tCurrent time
YSolution vector
YdotTime-derivative of solution vector
aShift
AJacobian matrix
BPreconditioning matrix
ctxtApplication context

Definition at line 52 of file PetscIJacobian.cpp.

60 {
61  PETScContext* context = (PETScContext*) ctxt;
62 
63  PetscFunctionBegin;
64  for (int ns = 0; ns < context->nsims; ns++) {
65  ((SimulationObject*)context->simobj)[ns].solver.count_IJacobian++;
66  }
67  context->shift = a;
68  context->waqt = t;
69  /* Construct preconditioning matrix */
70  if (context->flag_use_precon) PetscComputePreconMatImpl(B,Y,context);
71 
72  PetscFunctionReturn(0);
73 }
int PetscComputePreconMatImpl(Mat, Vec, void *)
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
PetscErrorCode PetscJacobianFunction_JFNK ( Mat  Jacobian,
Vec  Y,
Vec  F 
)

Computes the action of the Jacobian on a vector: See documentation for PetscIJacobian() 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 time integration of the ODE given by

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

we have that

\begin{equation} \mathcal{F}\left({\bf U}\right) \equiv \dot{\bf U} - {\bf F}\left({\bf U}\right) \Rightarrow {\bf J}\left({\bf U}\right) = \left[\alpha{\bf I} - \frac {\partial {\bf F}\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 F}\left({\bf U}_0+\epsilon{\bf y}\right)-{\bf F}\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.

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.
  • 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
JacobianJacobian matrix
YInput vector
FOutput vector (Jacobian times input vector)

Definition at line 111 of file PetscIJacobian.cpp.

115 {
116  PETScContext* context(nullptr);
117 
118  PetscFunctionBegin;
119 
120  MatShellGetContext(Jacobian,&context);
121  SimulationObject* sim = (SimulationObject*) context->simobj;
122  int nsims = context->nsims;
123 
124  double normY;
125  VecNorm(Y,NORM_2,&normY);
126 
127  if (normY < 1e-16) {
128 
129  /* F = 0 */
130  VecZeroEntries(F);
131  /* [J]Y = aY - F(Y) */
132  VecAXPBY(F,context->shift,0,Y);
133 
134  } else {
135 
136  double epsilon = context->jfnk_eps / normY;
137  double t = context->waqt; /* current stage/step time */
138 
139  for (int ns = 0; ns < nsims; ns++) {
140 
141  HyPar* solver = &(sim[ns].solver);
142  MPIVariables* mpi = &(sim[ns].mpi);
143  solver->count_IJacFunction++;
144 
145  int size = solver->npoints_local_wghosts;
146 
147  double *u = solver->u;
148  double *uref = solver->uref;
149  double *rhsref = solver->rhsref;
150  double *rhs = solver->rhs;
151 
152  /* copy solution from PETSc vector */
153  TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
154  _ArrayAYPX_(uref,epsilon,u,size*solver->nvars);
155  /* apply boundary conditions and exchange data over MPI interfaces */
156  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
158  solver->nvars,
159  solver->dim_local,
160  solver->ghosts,
161  mpi,
162  u );
163 
164  /* Evaluate hyperbolic, parabolic and source terms and the RHS for U+dU */
165  _ArraySetValue_(rhs,size*solver->nvars,0.0);
166  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
167  solver->FFunction,solver->Upwind);
168  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
169  solver->ParabolicFunction (solver->par,u,solver,mpi,t);
170  _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
171  solver->SourceFunction (solver->source,u,solver,mpi,t);
172  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
173 
174  _ArrayAXPY_(rhsref,-1.0,rhs,size*solver->nvars);
175  /* Transfer RHS to PETSc vector */
176  TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
177  }
178 
179  /* [J]Y = aY - F(Y) */
180  VecAXPBY(F,context->shift,(-1.0/epsilon),Y);
181 
182  }
183 
184  PetscFunctionReturn(0);
185 }
double * hyp
Definition: hypar.h:119
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int * dim_local
Definition: hypar.h:37
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
double * par
Definition: hypar.h:122
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
double * rhsref
Definition: hypar.h:398
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
double * uref
Definition: hypar.h:397
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int count_IJacFunction
Definition: hypar.h:422
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
#define _ArrayAYPX_(x, a, y, size)
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399
PetscErrorCode PetscJacobianFunction_Linear ( Mat  Jacobian,
Vec  Y,
Vec  F 
)

Computes the action of the Jacobian on a vector: See documentation for PetscIJacobian() 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 time integration of a linear ODE given by

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

we have that

\begin{equation} \mathcal{F}\left({\bf U}\right) \equiv \dot{\bf U} - {\bf F}\left({\bf U}\right) \Rightarrow {\bf J}\left({\bf U}\right) = \left[\alpha{\bf I} - \frac {\partial {\bf F}\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 F}\left({\bf U}_0+{\bf y}\right)-{\bf F}\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.

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
JacobianJacobian matrix
YInput vector
FOutput vector (Jacobian times input vector

Definition at line 236 of file PetscIJacobian.cpp.

240 {
241  PETScContext* context(nullptr);
242 
243  PetscFunctionBegin;
244 
245  MatShellGetContext(Jacobian,&context);
246  SimulationObject* sim = (SimulationObject*) context->simobj;
247  int nsims = context->nsims;
248 
249  double normY;
250  VecNorm(Y,NORM_2,&normY);
251 
252  if (normY < 1e-16) {
253 
254  /* F = 0 */
255  VecZeroEntries(F);
256  /* [J]Y = aY - F(Y) */
257  VecAXPBY(F,context->shift,0,Y);
258 
259  } else {
260 
261  double t = context->waqt; /* current stage/step time */
262 
263  for (int ns = 0; ns < nsims; ns++) {
264 
265  HyPar* solver = &(sim[ns].solver);
266  MPIVariables* mpi = &(sim[ns].mpi);
267  solver->count_IJacFunction++;
268 
269  int size = solver->npoints_local_wghosts;
270 
271  double *u = solver->u;
272  double *uref = solver->uref;
273  double *rhsref = solver->rhsref;
274  double *rhs = solver->rhs;
275 
276  /* copy solution from PETSc vector */
277  TransferVecFromPETSc(u,Y,context,ns,context->offsets[ns]);
278  _ArrayAYPX_(uref,1.0,u,size*solver->nvars);
279  /* apply boundary conditions and exchange data over MPI interfaces */
280  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,t);
282  solver->nvars,
283  solver->dim_local,
284  solver->ghosts,
285  mpi,
286  u );
287 
288  /* Evaluate hyperbolic, parabolic and source terms and the RHS for U+dU */
289  _ArraySetValue_(rhs,size*solver->nvars,0.0);
290  solver->HyperbolicFunction( solver->hyp,u,solver,mpi,t,0,
291  solver->FFunction,solver->Upwind);
292  _ArrayAXPY_(solver->hyp,-1.0,rhs,size*solver->nvars);
293  solver->ParabolicFunction (solver->par,u,solver,mpi,t);
294  _ArrayAXPY_(solver->par, 1.0,rhs,size*solver->nvars);
295  solver->SourceFunction (solver->source,u,solver,mpi,t);
296  _ArrayAXPY_(solver->source, 1.0,rhs,size*solver->nvars);
297 
298  _ArrayAXPY_(rhsref,-1.0,rhs,size*solver->nvars);
299  /* Transfer RHS to PETSc vector */
300  TransferVecToPETSc(rhs,F,context,ns,context->offsets[ns]);
301  }
302 
303  /* [J]Y = aY - F(Y) */
304  VecAXPBY(F,context->shift,-1.0,Y);
305 
306  }
307 
308  PetscFunctionReturn(0);
309 }
double * hyp
Definition: hypar.h:119
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int * dim_local
Definition: hypar.h:37
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
int ghosts
Definition: hypar.h:52
double * par
Definition: hypar.h:122
double * source
Definition: hypar.h:125
int(* SourceFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:259
double * rhsref
Definition: hypar.h:398
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
double * uref
Definition: hypar.h:397
int(* HyperbolicFunction)(double *, double *, void *, void *, double, int, int(*)(double *, double *, int, void *, double), int(*)(double *, double *, double *, double *, double *, double *, int, void *, double))
Definition: hypar.h:250
int count_IJacFunction
Definition: hypar.h:422
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
#define _ArrayAYPX_(x, a, y, size)
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399
int PetscGlobalDOF ( void *  c)

Compute the global DOF index for all the grid points: The "global DOF index" is the component number (or block component number for HyPar::nvars > 1) of a grid point in the global solution vector. It is also the row number (or block row number) of the grid point in the global matrix representing, for example, the Jacobian of the right-hand-side.

PETScContext::globalDOF is an integer array with the same layout as the solution array HyPar::u (but with one component) containing the global DOF index for the corresponding grid points. It has the same number of ghost points as HyPar::u.

  • This array is initialized to -1.
  • The global DOF indices are computed for all non-ghost grid points.
  • If any boundaries are periodic, periodic boundary conditions are applied to fill the appropriate ghost points.
  • Ghost points corresponding to internal (MPI) boundaries are filled using MPIExchangeBoundariesnD().
  • Thus, ghost points corresponding to physical, non-periodic boundaries retain the initial value of -1.
Parameters
cObject of type PETScContext

Definition at line 69 of file PetscGlobalDOF.cpp.

70 {
71  PETScContext* ctxt = (PETScContext*) c;
73  int nsims = ctxt->nsims;
74 
75  ctxt->globalDOF.resize(nsims, nullptr);
76 
77  /* compute MPI offset */
78  std::vector<int> local_sizes(ctxt->nproc ,0);
79  local_sizes[ctxt->rank] = ctxt->npoints;
80  MPIMax_integer(local_sizes.data(),local_sizes.data(),ctxt->nproc,&sim[0].mpi.world);
81 
82  int MPIOffset = 0;
83  for (int i=0; i<ctxt->rank; i++) MPIOffset += local_sizes[i];
84 
85  int simOffset = 0;
86  for (int ns = 0; ns < nsims; ns++) {
87 
88  HyPar* solver = &(sim[ns].solver);
89  MPIVariables* mpi = &(sim[ns].mpi);
90 
91  int *dim = solver->dim_local,
92  ndims = solver->ndims,
93  ghosts = solver->ghosts,
94  npoints = solver->npoints_local,
95  npoints_wg = solver->npoints_local_wghosts,
96  nv = ndims + 1, i;
97 
98  ctxt->globalDOF[ns] = (double*) calloc( npoints_wg, sizeof(double) );
99  _ArraySetValue_(ctxt->globalDOF[ns], npoints_wg, -1.0);
100 
101  for (int i = 0; i < npoints; i++) {
102  int p = (ctxt->points[ns]+i*nv)[ndims];
103  ctxt->globalDOF[ns][p] = (double) (i + simOffset + MPIOffset);
104  }
105 
106  for (int i=0; i<ndims; i++) {
107  if (solver->isPeriodic[i]) {
108  ApplyPeriodicity(i,ndims,dim,ghosts,ctxt->globalDOF[ns]);
109  }
110  }
111  MPIExchangeBoundariesnD(ndims,1,dim,ghosts,mpi,ctxt->globalDOF[ns]);
112 
113  simOffset += npoints;
114  }
115 
116  return 0;
117 }
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
int npoints_local
Definition: hypar.h:42
int MPIMax_integer(int *, int *, int, void *)
Definition: MPIMax.c:15
int * dim_local
Definition: hypar.h:37
std::vector< int * > points
static int ApplyPeriodicity(int dir, int ndims, int *size, int ghosts, double *phi)
int ghosts
Definition: hypar.h:52
MPI_Comm world
int * isPeriodic
Definition: hypar.h:162
std::vector< double * > globalDOF
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int PetscCleanup ( void *  obj)

Clean up allocations in the PETSc interface

Parameters
objObject of type PETScContext

Definition at line 12 of file PetscCleanup.cpp.

13 {
14  PETScContext *ctxt = (PETScContext*) obj;
15  for (int i = 0; i < ctxt->globalDOF.size(); i++) {
16  free(ctxt->globalDOF[i]);
17  }
18  ctxt->globalDOF.clear();
19  for (int i = 0; i < ctxt->points.size(); i++) {
20  free(ctxt->points[i]);
21  }
22  ctxt->points.clear();
23  if (ctxt->offsets) free(ctxt->offsets);
24  return(0);
25 }
std::vector< int * > points
std::vector< double * > globalDOF
Structure containing the variables for time-integration with PETSc.
int PetscCreatePointList ( void *  obj)

Create a list of computational points for each simulation domain: This is a list of all the grid points on which the PDE is solved. Thus, it is the total number of grid points minus the ghost points and blanked out points.

Note: this list is local, not global.

Parameters
objObject of type PETScContext

Definition at line 25 of file PetscCreatePointList.cpp.

26 {
27  PETScContext* ctxt = (PETScContext*) obj;
29 
30  int nsims = ctxt->nsims;
31  int ndims = sim[0].solver.ndims;
32 
33  /* count the number of computational points */
34  ctxt->npoints = 0;
35  ctxt->ndofs = 0;
36  ctxt->offsets = (int*) calloc(nsims, sizeof(int));
37  for (int ns = 0; ns < nsims; ns++) {
38  ctxt->npoints += sim[ns].solver.npoints_local;
39  ctxt->offsets[ns] = ctxt->ndofs;
40  ctxt->ndofs += (sim[ns].solver.npoints_local*sim[ns].solver.nvars);
41  }
42 
43  int nv = ndims+1;
44  ctxt->points.resize(nsims, nullptr);
45 
46  for (int ns = 0; ns < nsims; ns++) {
47 
48  int npoints = sim[ns].solver.npoints_local;
49  ctxt->points[ns] = (int*) calloc (npoints*nv,sizeof(int));
50 
51  const int* const dim( sim[ns].solver.dim_local );
52  const int ghosts( sim[ns].solver.ghosts );
53 
54  int done = 0, i = 0;
55  std::vector<int> index(ndims,0);
56  while (!done) {
57  int p; _ArrayIndex1D_(ndims, dim, index.data(), ghosts, p);
58  _ArrayCopy1D_(index.data(), (ctxt->points[ns]+i*nv), ndims);
59  (ctxt->points[ns]+i*nv)[ndims] = p;
60  _ArrayIncrementIndex_(ndims,dim,index,done);
61  i++;
62  }
63 
64  if (i != npoints) {
65  fprintf(stderr,"Error in PetscCreatePointList() on rank %d:\n", sim[ns].mpi.rank);
66  fprintf(stderr,"Inconsistency in point count - %d, %d.\n",
67  i, npoints);
68  return 1;
69  }
70  }
71 
72  int global_npoints;
73  int global_ndofs;
74  MPISum_integer(&global_npoints,&(ctxt->npoints),1,&sim[0].mpi.world);
75  MPISum_integer(&global_ndofs,&(ctxt->ndofs),1,&sim[0].mpi.world);
76  if (!ctxt->rank) {
77  printf("PETSc: total number of computational points is %d.\n",global_npoints);
78  printf("PETSc: total number of computational DOFs is %d.\n",global_ndofs);
79  }
80 
81  return 0;
82 }
int MPISum_integer(int *, int *, int, void *)
Definition: MPISum.c:16
int npoints_local
Definition: hypar.h:42
#define _ArrayIncrementIndex_(N, imax, i, done)
std::vector< int * > points
#define _ArrayIndex1D_(N, imax, i, ghost, index)
MPI_Comm world
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int ndims
Definition: hypar.h:26
PetscErrorCode PetscPreStage ( TS  ts,
PetscReal  waqt 
)

Function called before a stage in multi-stage time-integration methods

Parameters
tsTime integration object
waqtCurrent simulation time

Definition at line 14 of file PetscPreStage.cpp.

16 {
17  PetscFunctionBegin;
18  PetscFunctionReturn(0);
19 }
PetscErrorCode PetscPostStage ( TS  ts,
PetscReal  stagetime,
PetscInt  stageindex,
Vec *  Y 
)

Function called after every stage in a multi-stage time-integration method

Parameters
tsTime integrator of PETSc type TS
stagetimeCurrent stage time
stageindexStage
YStage solutions (all stages) - be careful what you access

Definition at line 16 of file PetscPostStage.cpp.

21 {
22  PETScContext* context(nullptr);
23 
24  PetscFunctionBegin;
25 
26  TSGetApplicationContext(ts,&context);
27  if (!context) {
28  fprintf(stderr,"Error in PetscPreTimeStep: Null context!\n");
29  return(1);
30  }
31  SimulationObject* sim = (SimulationObject*) context->simobj;
32  int nsims = context->nsims;
33 
34  TSType time_scheme;
35  TSGetType(ts,&time_scheme);
36 
37  TSGetTimeStep(ts,&(context->dt));
38  context->stage_index = stageindex;
39  if (context->stage_times.size() == stageindex) {
40  context->stage_times.push_back(stagetime/context->dt);
41  }
42 
43  for (int ns = 0; ns < nsims; ns++) {
44 
45  HyPar* solver = &(sim[ns].solver);
46  MPIVariables* mpi = &(sim[ns].mpi);
47 
48  /* get solution */
49  TransferVecFromPETSc(solver->u,Y[stageindex],context,ns,context->offsets[ns]);
50 
51  /* apply immersed boundaries */
52  solver->ApplyIBConditions(solver,mpi,solver->u,stagetime);
53 
54  /* If using a non-linear scheme with ARKIMEX methods,
55  compute the non-linear finite-difference operator */
56  if (!strcmp(time_scheme,TSARKIMEX)) {
57  solver->NonlinearInterp( solver->u,
58  solver,
59  mpi,
60  (double)stagetime,
61  solver->FFunction );
62  }
63 
64  /* Call any physics-specific post-stage function, if available */
65  if (solver->PostStage) {
66  solver->PostStage(solver->u,solver,mpi,stagetime);
67  }
68 
69  TransferVecToPETSc(solver->u,Y[stageindex],context,ns,context->offsets[ns]);
70 
71  }
72 
73  PetscFunctionReturn(0);
74 }
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* NonlinearInterp)(double *, void *, void *, double, int(*)(double *, double *, int, void *, double))
Definition: hypar.h:228
int(* PostStage)(double *, void *, void *, double)
Definition: hypar.h:336
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
Structure of MPI-related variables.
int(* ApplyIBConditions)(void *, void *, double *, double)
Definition: hypar.h:217
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
PetscErrorCode PetscPreTimeStep ( TS  ts)

Function called before a time step

Parameters
tsTime integration object

Definition at line 24 of file PetscPreTimeStep.cpp.

25 {
26  PETScContext* context(nullptr);
27 
28  PetscFunctionBegin;
29 
30  TSGetApplicationContext(ts,&context);
31  if (!context) {
32  fprintf(stderr,"Error in PetscPreTimeStep: Null context!\n");
33  return(1);
34  }
35  gettimeofday(&(context->iter_start_time),NULL);
36  SimulationObject* sim = (SimulationObject*) context->simobj;
37  int nsims = context->nsims;
38 
39  Vec Y;
40  TSGetSolution(ts,&Y);
41 
42  double waqt;
43  TSGetTime(ts,&waqt);
44  double dt;
45  TSGetTimeStep(ts,&dt);
46  int iter;
47  TSGetStepNumber(ts,&iter);
48 
49  context->dt = dt;
50  context->waqt = waqt;
51  context->t_start = waqt;
52 
53  TSType time_scheme;
54  TSGetType(ts,&time_scheme);
55 
56  for (int ns = 0; ns < nsims; ns++) {
57 
58  HyPar* solver = &(sim[ns].solver);
59  MPIVariables* mpi = &(sim[ns].mpi);
60 
61  /* get solution */
62  TransferVecFromPETSc(solver->u,Y,context,ns,context->offsets[ns]);
63 
64  /* save a copy of the solution to compute norm at end of time step */
65  _ArrayCopy1D_(solver->u,solver->u0,(solver->npoints_local_wghosts*solver->nvars));
66 
67  /* apply boundary conditions and exchange data over MPI interfaces */
68  solver->ApplyBoundaryConditions(solver,mpi,solver->u,NULL,waqt);
69  solver->ApplyIBConditions(solver,mpi,solver->u,waqt);
71  solver->nvars,
72  solver->dim_local,
73  solver->ghosts,
74  mpi,
75  solver->u );
76 
77  /* Call any physics-specific pre-step function */
78  if (solver->PreStep) solver->PreStep(solver->u,solver,mpi,waqt);
79 
80  /* If using a non-linear scheme with ARKIMEX methods,
81  compute the non-linear finite-difference operator */
82  if (!strcmp(time_scheme,TSARKIMEX)) {
83  solver->NonlinearInterp(solver->u,solver,mpi,waqt,solver->FFunction);
84  }
85 
86  /* set the step boundary flux integral value to zero */
87  _ArraySetValue_(solver->StepBoundaryIntegral,2*solver->ndims*solver->nvars,0.0);
88 
89  TransferVecToPETSc(solver->u,Y,context,ns,context->offsets[ns]);
90 
91  }
92 
93 
94  if (!iter) {
95  for (int ns = 0; ns < nsims; ns++) {
96  HyPar* solver = &(sim[ns].solver);
97  MPIVariables* mpi = &(sim[ns].mpi);
98  if (solver->PhysicsOutput) solver->PhysicsOutput(solver,mpi,waqt);
99  }
100  OutputSolution(sim, nsims,waqt);
101 #ifdef with_librom
102  context->op_times_arr.push_back(waqt);
103 #endif
104  }
105 
106 #ifdef with_librom
107  if ( (context->rom_mode == _ROM_MODE_TRAIN_)
108  && (iter%((libROMInterface*)context->rom_interface)->samplingFrequency() == 0) ) {
109  ((libROMInterface*)context->rom_interface)->takeSample( sim, waqt );
110  }
111 #endif
112 
113  PetscFunctionReturn(0);
114 }
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int * dim_local
Definition: hypar.h:37
double * u0
Definition: hypar.h:396
int(* PhysicsOutput)(void *, void *, double)
Definition: hypar.h:347
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int ghosts
Definition: hypar.h:52
int(* NonlinearInterp)(double *, void *, void *, double, int(*)(double *, double *, int, void *, double))
Definition: hypar.h:228
double * StepBoundaryIntegral
Definition: hypar.h:384
Class implementing interface with libROM.
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int OutputSolution(void *, int, double)
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int(* ApplyBoundaryConditions)(void *, void *, double *, double *, double)
Definition: hypar.h:214
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
int(* PreStep)(double *, void *, void *, double)
Definition: hypar.h:339
Structure of MPI-related variables.
#define _ROM_MODE_TRAIN_
int(* ApplyIBConditions)(void *, void *, double *, double)
Definition: hypar.h:217
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
PetscErrorCode PetscPostTimeStep ( TS  ts)

Function called after every time step

Parameters
tsTime integrator object

Definition at line 21 of file PetscPostTimeStep.cpp.

22 {
23  PETScContext* context(nullptr);
24 
25  PetscFunctionBegin;
26 
27  TSGetApplicationContext(ts,&context);
28  if (!context) {
29  fprintf(stderr,"Error in PetscPreTimeStep: Null context!\n");
30  return(1);
31  }
32  SimulationObject* sim = (SimulationObject*) context->simobj;
33  int nsims = context->nsims;
34 
35  Vec Y;
36  TSGetSolution(ts,&Y);
37 
38  double waqt;
39  TSGetTime(ts,&waqt);
40  int iter;
41  TSGetStepNumber(ts,&iter);
42  double dt;
43  TSGetTimeStep(ts,&dt);
44 
45  context->tic++;
46 
47  double max_cfl = 0.0;
48  double total_norm = 0.0;
49  for (int ns = 0; ns < nsims; ns++) {
50 
51  HyPar* solver = &(sim[ns].solver);
52  MPIVariables* mpi = &(sim[ns].mpi);
53 
54  /* get solution */
55  TransferVecFromPETSc(solver->u,Y,context,ns,context->offsets[ns]);
56 
57  /* Call any physics-specific post-step function */
58  if (solver->PostStep) solver->PostStep(solver->u,solver,mpi,waqt,iter);
59 
60  /* Calculate CFL and diffusion number */
61  double local_max_cfl = -1.0, global_max_cfl = -1.0;
62  if (solver->ComputeCFL) local_max_cfl = solver->ComputeCFL(solver,mpi,dt,waqt);
63  MPIMax_double(&global_max_cfl ,&local_max_cfl ,1,&mpi->world);
64  if (global_max_cfl > max_cfl) max_cfl = global_max_cfl;
65 
66  /* Calculate norm of the change in the solution for this time step */
67  _ArrayAXPY_(solver->u,-1.0,solver->u0,(solver->npoints_local_wghosts*solver->nvars));
68  double sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
69  solver->ghosts,solver->index,solver->u0);
70  double global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
71  double norm = sqrt((global_sum/(double)solver->npoints_global));
72  total_norm += norm;
73 
74  if (!strcmp(solver->ConservationCheck,"yes")) {
75  /* calculate volume integral of the solution at this time step */
76  solver->VolumeIntegralFunction(solver->VolumeIntegral,solver->u,solver,mpi);
77  /* calculate surface integral of the flux at this time step */
78  solver->BoundaryIntegralFunction(solver,mpi);
79  /* calculate the conservation error at this time step */
80  solver->CalculateConservationError(solver,mpi);
81  }
82 
83  }
84 
85  gettimeofday(&context->iter_end_time,NULL);
86  long long walltime;
87  walltime = ( (context->iter_end_time.tv_sec * 1000000 + context->iter_end_time.tv_usec)
88  - (context->iter_start_time.tv_sec * 1000000 + context->iter_start_time.tv_usec));
89  context->iter_wctime = (double) walltime / 1000000.0;
90  context->ti_runtime += context->iter_wctime;
91 
92  if ((!context->rank) && (iter%sim[0].solver.screen_op_iter == 0)) {
93 
94  if (nsims > 1) {
95  printf("--\n");
96  printf("iter=%7d, t=%1.3e\n", iter, waqt);
97  printf(" dt=%1.3E ", dt);
98  printf(" CFL=%1.3E, ", max_cfl);
99  printf(" norm=%1.4E\n", total_norm);
100  printf(" wctime=%1.1E (s)\n",context->iter_wctime);
101  } else {
102  printf("iter=%7d ", iter);
103  printf("dt=%1.3E ", dt);
104  printf("t=%1.3E ", waqt);
105  printf("CFL=%1.3E ", max_cfl );
106  printf("norm=%1.4E ", total_norm);
107  printf("wctime: %1.1E (s) ",context->iter_wctime);
108  }
109 
110  /* calculate and print conservation error */
111  if (!strcmp(sim[0].solver.ConservationCheck,"yes")) {
112 
113  double error = 0;
114  for (int ns = 0; ns < nsims; ns++) {
115  for (int v=0; v<sim[ns].solver.nvars; v++) {
116  error += (sim[ns].solver.ConservationError[v] * sim[ns].solver.ConservationError[v]);
117  }
118  }
119  error = sqrt(error);
120  printf(" cons_err=%1.4E\n", error);
121 
122  } else {
123 
124  if (nsims == 1) printf("\n");
125 
126  }
127 
128  /* print physics-specific info, if available */
129  for( int ns = 0; ns < nsims; ns++) {
130  if (sim[ns].solver.PrintStep) {
131  if (nsims > 1) printf("Physics-specific output for domain %d:\n", ns);
132  sim[ns].solver.PrintStep( &(sim[ns].solver),
133  &(sim[ns].mpi),
134  waqt );
135  }
136  }
137 
138  if (nsims > 1) {
139  printf("--\n");
140  printf("\n");
141  }
142 
143  }
144 
145  /* Write intermediate solution to file */
146  if (iter%sim[0].solver.file_op_iter == 0) {
147  for (int ns = 0; ns < nsims; ns++) {
148  HyPar* solver = &(sim[ns].solver);
149  MPIVariables* mpi = &(sim[ns].mpi);
150  if (solver->PhysicsOutput) solver->PhysicsOutput(solver,mpi,waqt);
151  }
152  OutputSolution(sim,nsims,waqt);
153 #ifdef with_librom
154  context->op_times_arr.push_back(waqt);
155 #endif
156  context->tic=0;
157  }
158 
159  PetscFunctionReturn(0);
160 }
int npoints_local_wghosts
Definition: hypar.h:42
double(* ComputeCFL)(void *, void *, double, double)
Definition: hypar.h:269
double * u
Definition: hypar.h:116
int npoints_global
Definition: hypar.h:40
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
int * dim_local
Definition: hypar.h:37
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
double * u0
Definition: hypar.h:396
int(* PhysicsOutput)(void *, void *, double)
Definition: hypar.h:347
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
int ghosts
Definition: hypar.h:52
int screen_op_iter
Definition: hypar.h:168
MPI_Comm world
char ConservationCheck[_MAX_STRING_SIZE_]
Definition: hypar.h:376
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
double * VolumeIntegral
Definition: hypar.h:378
int(* BoundaryIntegralFunction)(void *, void *)
Definition: hypar.h:390
int(* VolumeIntegralFunction)(double *, double *, void *, void *)
Definition: hypar.h:388
int nvars
Definition: hypar.h:29
int file_op_iter
Definition: hypar.h:171
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
int OutputSolution(void *, int, double)
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int(* PrintStep)(void *, void *, double)
Definition: hypar.h:344
int ndims
Definition: hypar.h:26
int * index
Definition: hypar.h:102
double * ConservationError
Definition: hypar.h:374
int(* CalculateConservationError)(void *, void *)
Definition: hypar.h:392
Structure of MPI-related variables.
#define _ArrayAXPY_(x, a, y, size)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int(* PostStep)(double *, void *, void *, double, int)
Definition: hypar.h:341
int PetscComputePreconMatIMEX ( Mat  Pmat,
Vec  Y,
void *  ctxt 
)

Compute and assemble the preconditioning matrix for the implicit-explicit (IMEX) time integration of the governing equations: Right now, it just calls PetscComputePreconMatImpl()

Parameters
PmatPreconditioning matrix to construct
YSolution vector
ctxtApplication context

Definition at line 17 of file PetscComputePreconMatIMEX.cpp.

20 {
21  /* Same implementation as PetscComputePreconMatImpl() */
22  return(PetscComputePreconMatImpl(Pmat,Y,ctxt));
23 }
int PetscComputePreconMatImpl(Mat, Vec, void *)
int PetscComputePreconMatImpl ( Mat  Pmat,
Vec  Y,
void *  ctxt 
)

Compute and assemble the preconditioning matrix for the implicit 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 L}\left({\bf U}\right) \Rightarrow \frac {d{\bf U}}{dt} - {\bf L}\left({\bf U}\right) = 0, \end{equation}

where \({\bf L}\) is the spatially discretized right-hand-side, and \({\bf U}\) represents the entire solution vector.

The Jacobian is thus given by:

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

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

Let \(\mathcal{D}_{\rm hyp}\) and \(\mathcal{D}_{\rm par}\) represent the spatial discretization methods for the hyperbolic and parabolic terms. Thus, the Jacobian can be written as follows:

\begin{equation} {\bf J} = \left[\alpha{\bf I} - \left(\mathcal{D}_{\rm hyp}\left\{\frac {\partial {\bf F}_{\rm hyp}} {\partial {\bf u}}\right\} + \mathcal{D}_{\rm par}\left\{\frac {\partial {\bf F}_{\rm par}} {\partial {\bf u}}\right\}\right)\right] \end{equation}

The preconditioning matrix is usually a close approximation of the actual Jacobian matrix, where the actual Jacobian may be too expensive to evaluate and assemble. In this function, the preconditioning matrix is the following approximation of the actual Jacobian:

\begin{equation} {\bf J}_p = \left[\alpha{\bf I} - \left(\mathcal{D}^{\left(l\right)}_{\rm hyp}\left\{\frac {\partial {\bf F}_{\rm hyp}} {\partial {\bf u}}\right\} + \mathcal{D}^{\left(l\right)}_{\rm par}\left\{\frac {\partial {\bf F}_{\rm par}} {\partial {\bf u}}\right\}\right) \right] \approx {\bf J}, \end{equation}

where \(\mathcal{D}^{\left(l\right)}_{\rm hyp,par}\) represent lower order discretizations of the hyperbolic and parabolic terms. The matrix \({\bf J}_p\) is provided to the preconditioner.

Note that the specific physics model provides the following functions:

  • HyPar::JFunction computes \(\partial {\bf F}_{\rm hyp}/ \partial {\bf u}\) at a grid point.
  • HyPar::KFunction computes \(\partial {\bf F}_{\rm par}/ \partial {\bf u}\) at a grid point.

Currently, this function doesn't include the source term.

Parameters
PmatPreconditioning matrix to construct
YSolution vector
ctxtApplication context

Definition at line 59 of file PetscComputePreconMatImpl.cpp.

62 {
63  PetscErrorCode ierr;
64  PETScContext* context = (PETScContext*) ctxt;
65  SimulationObject* sim = (SimulationObject*) context->simobj;
66  int nsims = context->nsims;
67 
68  PetscFunctionBegin;
69  /* initialize preconditioning matrix to zero */
70  MatZeroEntries(Pmat);
71 
72  /* copy solution from PETSc vector */
73  for (int ns = 0; ns < nsims; ns++) {
74 
75  TransferVecFromPETSc( sim[ns].solver.u,
76  Y,
77  context,
78  ns,
79  context->offsets[ns]);
80 
81  HyPar* solver( &(sim[ns].solver) );
82  MPIVariables* mpi( &(sim[ns].mpi) );
83 
84  int ndims = solver->ndims,
85  nvars = solver->nvars,
86  npoints = solver->npoints_local,
87  ghosts = solver->ghosts,
88  *dim = solver->dim_local,
89  *isPeriodic = solver->isPeriodic,
90  *points = context->points[ns],
91  index[ndims],indexL[ndims],indexR[ndims],
92  rows[nvars],cols[nvars];
93 
94  double *u = solver->u,
95  *iblank = solver->iblank,
96  dxinv, values[nvars*nvars];
97 
98  /* apply boundary conditions and exchange data over MPI interfaces */
99  solver->ApplyBoundaryConditions(solver,mpi,u,NULL,context->waqt);
100  MPIExchangeBoundariesnD(ndims,nvars,dim,ghosts,mpi,u);
101 
102  /* loop through all computational points */
103  for (int n = 0; n < npoints; n++) {
104  int *this_point = points + n*(ndims+1);
105  int p = this_point[ndims];
106  int index[ndims]; _ArrayCopy1D_(this_point,index,ndims);
107 
108  double iblank = solver->iblank[p];
109 
110  /* compute the contributions from the hyperbolic flux derivatives along each dimension */
111  if (solver->JFunction) {
112  for (int dir = 0; dir < ndims; dir++) {
113 
114  /* compute indices and global 1D indices for left and right neighbors */
115  _ArrayCopy1D_(index,indexL,ndims); indexL[dir]--;
116  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
117 
118  _ArrayCopy1D_(index,indexR,ndims); indexR[dir]++;
119  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
120 
121  int pg, pgL, pgR;
122  pg = (int) context->globalDOF[ns][p];
123  pgL = (int) context->globalDOF[ns][pL];
124  pgR = (int) context->globalDOF[ns][pR];
125 
126  /* Retrieve 1/delta-x at this grid point */
127  _GetCoordinate_(dir,index[dir],dim,ghosts,solver->dxinv,dxinv);
128 
129  /* diagonal element */
130  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pg + v; }
131  solver->JFunction(values,(u+nvars*p),solver->physics,dir,nvars,0);
132  _ArrayScale1D_(values,(dxinv*iblank),(nvars*nvars));
133  MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
134 
135  /* left neighbor */
136  if (pgL >= 0) {
137  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgL + v; }
138  solver->JFunction(values,(u+nvars*pL),solver->physics,dir,nvars,1);
139  _ArrayScale1D_(values,(-dxinv*iblank),(nvars*nvars));
140  MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
141  }
142 
143  /* right neighbor */
144  if (pgR >= 0) {
145  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgR + v; }
146  solver->JFunction(values,(u+nvars*pR),solver->physics,dir,nvars,-1);
147  _ArrayScale1D_(values,(-dxinv*iblank),(nvars*nvars));
148  MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
149  }
150  }
151  }
152 
153  /* compute the contributions from the parabolic term derivatives along each dimension */
154  if (solver->KFunction) {
155  for (int dir = 0; dir < ndims; dir++) {
156 
157  /* compute indices and global 1D indices for left and right neighbors */
158  _ArrayCopy1D_(index,indexL,ndims); indexL[dir]--;
159  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
160 
161  _ArrayCopy1D_(index,indexR,ndims); indexR[dir]++;
162  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
163 
164  int pg, pgL, pgR;
165  pg = (int) context->globalDOF[ns][p];
166  pgL = (int) context->globalDOF[ns][pL];
167  pgR = (int) context->globalDOF[ns][pR];
168 
169  /* Retrieve 1/delta-x at this grid point */
170  _GetCoordinate_(dir,index[dir],dim,ghosts,solver->dxinv,dxinv);
171 
172  /* diagonal element */
173  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pg + v; }
174  solver->KFunction(values,(u+nvars*p),solver->physics,dir,nvars);
175  _ArrayScale1D_(values,(-2*dxinv*dxinv*iblank),(nvars*nvars));
176  MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
177 
178  /* left neighbor */
179  if (pgL >= 0) {
180  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgL + v; }
181  solver->KFunction(values,(u+nvars*pL),solver->physics,dir,nvars);
182  _ArrayScale1D_(values,(dxinv*dxinv*iblank),(nvars*nvars));
183  MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
184  }
185 
186  /* right neighbor */
187  if (pgR >= 0) {
188  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgR + v; }
189  solver->KFunction(values,(u+nvars*pR),solver->physics,dir,nvars);
190  _ArrayScale1D_(values,(dxinv*dxinv*iblank),(nvars*nvars));
191  MatSetValues(Pmat,nvars,rows,nvars,cols,values,ADD_VALUES);
192  }
193  }
194  }
195  }
196  }
197 
198  MatAssemblyBegin(Pmat,MAT_FINAL_ASSEMBLY);
199  MatAssemblyEnd (Pmat,MAT_FINAL_ASSEMBLY);
200 
201  MatShift(Pmat,context->shift);
202  PetscFunctionReturn(0);
203 }
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
std::vector< int * > points
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
#define _ArrayScale1D_(x, a, size)
#define _ArrayIndex1D_(N, imax, i, ghost, index)
std::vector< double * > globalDOF
#define _ArrayCopy1D_(x, y, size)
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int PetscJacobianMatNonzeroEntriesImpl ( Mat  Amat,
int  width,
void *  ctxt 
)

Set non-zero entries of the Jacobian matrix for the implicit 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 L}\left({\bf U}\right) \Rightarrow \frac {d{\bf U}}{dt} - {\bf L}\left({\bf U}\right) = 0, \end{equation}

where \({\bf L}\) is the spatially discretized right-hand-side, and \({\bf U}\) represents the entire solution vector.

The Jacobian is thus given by:

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

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

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
AmatMatrix
widthStencil width
ctxtApplication context

Definition at line 38 of file PetscJacobianMatNonzeroEntriesImpl.cpp.

41 {
42  PETScContext* context = (PETScContext*) ctxt;
43  SimulationObject* sim = (SimulationObject*) context->simobj;
44 
45  PetscFunctionBegin;
46  int nsims = context->nsims;
47  /* initialize matrix to zero */
48  MatZeroEntries(Amat);
49 
50  for (int ns = 0; ns < nsims; ns++) {
51 
52  HyPar* solver( &(sim[ns].solver) );
53  MPIVariables* mpi( &(sim[ns].mpi) );
54 
55  int ndims = solver->ndims,
56  nvars = solver->nvars,
57  npoints = solver->npoints_local,
58  ghosts = solver->ghosts,
59  *dim = solver->dim_local,
60  *points = context->points[ns],
61  index[ndims],indexL[ndims],indexR[ndims],
62  rows[nvars],cols[nvars];
63 
64  std::vector<double> values(nvars*nvars, 1.0);
65 
66  /* loop through all computational points */
67  for (int n = 0; n < npoints; n++) {
68 
69  int *this_point = points + n*(ndims+1);
70  int p = this_point[ndims];
71  int index[ndims]; _ArrayCopy1D_(this_point,index,ndims);
72 
73  for (int dir = 0; dir < ndims; dir++) {
74 
75  int pg = (int) context->globalDOF[ns][p];
76  /* diagonal element */
77  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pg + v; }
78  MatSetValues(Amat,nvars,rows,nvars,cols,values.data(),ADD_VALUES);
79 
80  for (int d = 1; d <= width; d++) {
81 
82  /* left neighbor */
83  _ArrayCopy1D_(index,indexL,ndims);
84  indexL[dir] -= d;
85  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
86  int pgL = (int) context->globalDOF[ns][pL];
87  if (pgL >= 0) {
88  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgL + v; }
89  MatSetValues(Amat,nvars,rows,nvars,cols,values.data(),ADD_VALUES);
90  }
91 
92  _ArrayCopy1D_(index,indexR,ndims);
93  indexR[dir] += d;
94  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
95  int pgR = (int) context->globalDOF[ns][pR];
96  /* right neighbor */
97  if (pgR >= 0) {
98  for (int v=0; v<nvars; v++) { rows[v] = nvars*pg + v; cols[v] = nvars*pgR + v; }
99  MatSetValues(Amat,nvars,rows,nvars,cols,values.data(),ADD_VALUES);
100  }
101 
102  }
103 
104  }
105  }
106  }
107 
108  MatAssemblyBegin(Amat,MAT_FINAL_ASSEMBLY);
109  MatAssemblyEnd (Amat,MAT_FINAL_ASSEMBLY);
110 
111  PetscFunctionReturn(0);
112 }
std::vector< int * > points
#define _ArrayIndex1D_(N, imax, i, ghost, index)
std::vector< double * > globalDOF
#define _ArrayCopy1D_(x, y, size)
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
PetscErrorCode PetscTimeError ( TS  ts)

Function to compute any error estimates, if available

Compute the norms of the error estimate, if the PETSc time integrator has it (for example TSGLEE)

Parameters
tsTime integrator object of PETSc type TS

Definition at line 20 of file PetscError.cpp.

21 {
22  PetscErrorCode ierr;
23 
24  PETScContext* context(nullptr);
25  ierr = TSGetApplicationContext(ts,&context); CHKERRQ(ierr);
26  if (!context) {
27  fprintf(stderr,"Error in PetscError: Null context!\n");
28  return(1);
29  }
30 
31  SimulationObject* sim( (SimulationObject*)context->simobj );
32  int nsims( context->nsims );
33 
34  double dt;
35  ierr = TSGetTimeStep(ts,&dt); CHKERRQ(ierr);
36  TSType time_scheme;
37  ierr = TSGetType(ts,&time_scheme); CHKERRQ(ierr);
38  Vec Y;
39  ierr = TSGetSolution(ts,&Y); CHKERRQ(ierr);
40 
41  for (int ns = 0; ns < nsims; ns++) {
42  TransferVecFromPETSc( sim[ns].solver.u,
43  Y,
44  context,
45  ns,
46  context->offsets[ns] );
47  }
48 
49  if (std::string(time_scheme) == std::string(TSGLEE)) {
50 
51  Vec Z;
52  ierr = VecDuplicate(Y,&Z); CHKERRQ(ierr);
53  ierr = TSGetTimeError(ts,0,&Z);CHKERRQ(ierr);
54  for (int ns = 0; ns < nsims; ns++) {
55  TransferVecFromPETSc( sim[ns].solver.uref,
56  Z,
57  context,
58  ns,
59  context->offsets[ns] );
60  }
61  ierr = VecDestroy(&Z); CHKERRQ(ierr);
62 
63  for (int ns = 0; ns < nsims; ns++) {
64 
65  HyPar* solver = &(sim[ns].solver);
66  MPIVariables* mpi = &(sim[ns].mpi);
67 
68  int size = solver->npoints_local_wghosts * solver->nvars;
69  double sum = 0.0,
70  global_sum = 0.0,
71  *Uerr = solver->uref,
72  error[3] = {0,0,0};
73 
74  /* calculate solution norm for relative errors */
75  double sol_norm[3] = {0.0,0.0,0.0};
76  /* L1 */
77  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
78  solver->ghosts,solver->index,solver->u);
79  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
80  sol_norm[0] = global_sum/((double)solver->npoints_global);
81  /* L2 */
82  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
83  solver->ghosts,solver->index,solver->u);
84  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
85  sol_norm[1] = sqrt(global_sum/((double)solver->npoints_global));
86  /* Linf */
87  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
88  solver->ghosts,solver->index,solver->u);
89  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
90  sol_norm[2] = global_sum;
91 
92  /* calculate L1 norm of error */
93  sum = ArraySumAbsnD (solver->nvars,solver->ndims,solver->dim_local,
94  solver->ghosts,solver->index,Uerr);
95  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
96  error[0] = global_sum/((double)solver->npoints_global);
97  /* calculate L2 norm of error */
98  sum = ArraySumSquarenD(solver->nvars,solver->ndims,solver->dim_local,
99  solver->ghosts,solver->index,Uerr);
100  global_sum = 0; MPISum_double(&global_sum,&sum,1,&mpi->world);
101  error[1] = sqrt(global_sum/((double)solver->npoints_global));
102  /* calculate Linf norm of error */
103  sum = ArrayMaxnD (solver->nvars,solver->ndims,solver->dim_local,
104  solver->ghosts,solver->index,Uerr);
105  global_sum = 0; MPIMax_double(&global_sum,&sum,1,&mpi->world);
106  error[2] = global_sum;
107 
108  if ( (sol_norm[0] > _MACHINE_ZERO_)
109  && (sol_norm[1] > _MACHINE_ZERO_)
110  && (sol_norm[2] > _MACHINE_ZERO_) ) {
111  error[0] /= sol_norm[0];
112  error[1] /= sol_norm[1];
113  error[2] /= sol_norm[2];
114  }
115 
116  /* write to file */
117  if (!mpi->rank) {
118  std::string fname = "glm_err";
119  if (nsims > 1) {
120  char idx_string[10];
121  sprintf(idx_string, "%3d", ns);
122  fname += ("_" + std::string(idx_string));
123  }
124  fname += ".dat";
125 
126  FILE *out;
127  out = fopen(fname.c_str(),"w");
128  fprintf(out,"%1.16E %1.16E %1.16E %1.16E ",dt,error[0],error[1],error[2]);
129  fclose(out);
130  if (nsims > 1) {
131  printf( "Estimated time integration errors (GLM-GEE time-integration) for simulation domain %d:-\n",
132  ns);
133  } else {
134  printf("Estimated time integration errors (GLM-GEE time-integration):-\n");
135  }
136  printf(" L1 Error : %1.16E\n",error[0]);
137  printf(" L2 Error : %1.16E\n",error[1]);
138  printf(" Linfinity Error : %1.16E\n",error[2]);
139  }
140  }
141  }
142 
143  return 0;
144 }
int npoints_local_wghosts
Definition: hypar.h:42
double * u
Definition: hypar.h:116
int npoints_global
Definition: hypar.h:40
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
INLINE double ArraySumAbsnD(int, int, int *, int, int *, double *)
int * dim_local
Definition: hypar.h:37
int MPISum_double(double *, double *, int, void *)
Definition: MPISum.c:39
int MPIMax_double(double *, double *, int, void *)
Definition: MPIMax.c:38
int ghosts
Definition: hypar.h:52
MPI_Comm world
INLINE double ArraySumSquarenD(int, int, int *, int, int *, double *)
double * uref
Definition: hypar.h:397
int nvars
Definition: hypar.h:29
long sum(const std::vector< int > &a_iv)
Definition: std_vec_ops.h:18
Structure containing the variables for time-integration with PETSc.
#define _MACHINE_ZERO_
Definition: basic.h:26
Structure defining a simulation.
int ndims
Definition: hypar.h:26
int * index
Definition: hypar.h:102
Structure of MPI-related variables.
INLINE double ArrayMaxnD(int, int, int *, int, int *, double *)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
PetscErrorCode PetscSetInitialGuessROM ( SNES  snes,
Vec  X,
void *  ctxt 
)

Function to compute initial guess of nonlinear solve using libROM

Compute the initial guess for a nonlinear solve using a trained libROM reduced-order model.

Parameters
snesNonlinear solver object (see PETSc SNES)
XInitial guess vector
ctxtApplication context

Definition at line 18 of file PetscSetInitialGuessROM.cpp.

21 {
22  PETScContext* context = (PETScContext*) ctxt;
23  SimulationObject* sim = (SimulationObject*) context->simobj;
24  int nsims = context->nsims;
25 
26  PetscFunctionBegin;
27 
28  double stage_time = context->t_start;
29  if (context->stage_times.size() > (context->stage_index+1)) {
30  stage_time += context->stage_times[context->stage_index+1] * context->dt;
31  } else {
32  stage_time += context->dt;
33  }
34 
35  ((libROMInterface*)context->rom_interface)->predict(sim, stage_time);
36  if (!context->rank) {
37  printf( "libROM: Predicted ROM intial guess (time %1.4e), wallclock time: %f.\n",
38  stage_time, ((libROMInterface*)context->rom_interface)->predictWallclockTime() );
39  }
40 
41  for (int ns = 0; ns < nsims; ns++) {
42  TransferVecToPETSc(sim[ns].solver.u,X,context,ns,context->offsets[ns]);
43  }
44 
45 
46  PetscFunctionReturn(0);
47 }
std::vector< double > stage_times
Class implementing interface with libROM.
int TransferVecToPETSc(const double *const, Vec, void *, const int, const int)
Structure containing the variables for time-integration with PETSc.
Structure defining a simulation.