HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
PetscIJacobianIMEX.cpp
Go to the documentation of this file.
1 
6 #ifdef with_petsc
7 
8 #include <stdlib.h>
9 #include <basic.h>
10 #include <arrayfunctions.h>
11 #include <mpivars_cpp.h>
12 #include <simulation_object.h>
13 #include <petscinterface.h>
14 
15 #undef __FUNCT__
16 #define __FUNCT__ "PetscIJacobianIMEX"
17 
61 PetscErrorCode PetscIJacobianIMEX(TS ts,
62  PetscReal t,
63  Vec Y,
64  Vec Ydot,
65  PetscReal a,
66  Mat A,
67  Mat B,
68  void *ctxt )
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 }
83 
84 #undef __FUNCT__
85 #define __FUNCT__ "PetscJacobianFunctionIMEX_JFNK"
86 
128 PetscErrorCode PetscJacobianFunctionIMEX_JFNK(Mat Jacobian,
129  Vec Y,
130  Vec F )
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 }
235 
236 #undef __FUNCT__
237 #define __FUNCT__ "PetscJacobianFunctionIMEX_Linear"
238 
290 PetscErrorCode PetscJacobianFunctionIMEX_Linear(Mat Jacobian,
291  Vec Y,
292  Vec F )
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 }
397 
398 #endif
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
PetscErrorCode PetscIJacobianIMEX(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *)
double * u
Definition: hypar.h:116
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)
PetscErrorCode PetscJacobianFunctionIMEX_JFNK(Mat, Vec, Vec)
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
int PetscComputePreconMatIMEX(Mat, Vec, void *)
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
Simulation object.
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
Some basic definitions and macros.
int ndims
Definition: hypar.h:26
Contains macros and function definitions for common array operations.
PetscErrorCode PetscJacobianFunctionIMEX_Linear(Mat, Vec, Vec)
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)
C++ declarations for MPI-related functions.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double * rhs
Definition: hypar.h:399