HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
PetscComputePreconMatImpl.cpp File Reference

Contains the function to assemble the preconditioning matrix. More...

#include <stdio.h>
#include <arrayfunctions.h>
#include <simulation_object.h>
#include <mpivars_cpp.h>
#include <petscinterface.h>

Go to the source code of this file.

Macros

#define __FUNCT__   "PetscComputePreconMatImpl"
 

Functions

int PetscComputePreconMatImpl (Mat Pmat, Vec Y, void *ctxt)
 

Detailed Description

Contains the function to assemble the preconditioning matrix.

Author
Debojyoti Ghosh

Definition in file PetscComputePreconMatImpl.cpp.

Macro Definition Documentation

◆ __FUNCT__

#define __FUNCT__   "PetscComputePreconMatImpl"

Definition at line 15 of file PetscComputePreconMatImpl.cpp.

Function Documentation

◆ PetscComputePreconMatImpl()

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 }
std::vector< int * > points
Structure defining a simulation.
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
#define _ArrayIndex1D_(N, imax, i, ghost, index)
Structure containing the variables for time-integration with PETSc.
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
Structure of MPI-related variables.
#define _ArrayScale1D_(x, a, size)
std::vector< double * > globalDOF
#define _ArrayCopy1D_(x, y, size)
int TransferVecFromPETSc(double *const, const Vec, void *, const int, const int)