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

Initialize the linear advection-diffusion-reaction module. More...

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <bandedmatrix.h>
#include <physicalmodels/linearadr.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions

int LinearADRAdvectionField (void *, void *, int, int, int *)
 
double LinearADRComputeCFL (void *, void *, double, double)
 
double LinearADRComputeDiffNumber (void *, void *, double, double)
 
int LinearADRAdvection (double *, double *, int, void *, double)
 
int LinearADRDiffusionG (double *, double *, int, void *, double)
 
int LinearADRDiffusionH (double *, double *, int, int, void *, double)
 
int LinearADRReaction ()
 
int LinearADRUpwind (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int LinearADRCenteredFlux (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int LinearADRWriteAdvField (void *, void *, double)
 
int LinearADRAdvectionJacobian (double *, double *, void *, int, int, int)
 
int LinearADRDiffusionJacobian (double *, double *, void *, int, int)
 
int LinearADRInitialize (void *s, void *m)
 

Detailed Description

Initialize the linear advection-diffusion-reaction module.

Author
Debojyoti Ghosh

Definition in file LinearADRInitialize.c.

Function Documentation

int LinearADRAdvectionField ( void *  s,
void *  m,
int  idx,
int  nsims,
int *  dim_adv 
)

Set the advection field over the domain - reads the spatially-varying advection data from a file, if available. The array to store the field must already be allocated.

For a single simulation, it reads in the data from a file. The data must have the same grid dimensions as the solution.

For an ensemble simulation, it reads in the advection field from the file corresponding to the idx (index of this simulation). The data in this file must have the same grid dimensions as the domain it is being read in for.

For a sparse grids simulation, it reads in the advection field from the file with data that has the grid dimensions as the full grid. The field on the current grid is obtained by interpolation.

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables
idxIndex of this simulation
nsimsTotal number of simulations
dim_advGrid dimensions of the advection field stored in file

Definition at line 33 of file LinearADRAdvectionField.c.

39 {
40  HyPar *solver = (HyPar*) s;
41  MPIVariables *mpi = (MPIVariables*) m;
42  LinearADR *param = (LinearADR*) solver->physics;
43  double *adv = param->a;
44 
45  /* sanity checks */
46  if (param->constant_advection == 1) {
47  fprintf(stderr,"Error in LinearADRAdvectionField(): param->constant_advection == 1!\n");
48  return(1);
49  }
50  if (!strcmp(param->adv_filename,"none")) {
51  fprintf(stderr,"Error in LinearADRAdvectionField(): param->adv_filename is \"none\"!\n");
52  return(1);
53  }
54 
55  int d, done, flag;
56  int *dim = solver->dim_local;
57  int ghosts = solver->ghosts;
58 
59  /* number of advection field components at a grid point is
60  * number of spatial dimensions times number of solution
61  * components */
62  int adv_nvar = solver->ndims * solver->nvars;
63  if (param->adv_arr_size != adv_nvar*solver->npoints_local_wghosts) {
64  fprintf(stderr,"Error in LinearADRAdvectionField(): Incorrect adv_arr_size\n");
65  return(1);
66  }
67 
68  char fname_root[_MAX_STRING_SIZE_];
69  strcpy(fname_root, param->adv_filename);
70  if (idx >= 0) {
71  if (nsims > 1) {
72  char index[_MAX_STRING_SIZE_];
73  GetStringFromInteger(idx, index, (int)log10(nsims)+1);
74  strcat(fname_root, "_");
75  strcat(fname_root, index);
76  }
77  }
78 
79  /* read spatially-varying advection field from provided file */
80  if (dim_adv == NULL) {
81  int ierr = ReadArray( solver->ndims,
82  adv_nvar,
83  solver->dim_global,
84  solver->dim_local,
85  solver->ghosts,
86  solver,
87  mpi,
88  NULL,
89  adv,
90  fname_root,
91  &flag);
92  if (ierr) {
93  fprintf(stderr,"Error in LinearADRAdvectionField()\n");
94  fprintf(stderr," ReadArray() returned error!\n");
95  return ierr;
96  }
97  } else {
98  int ierr = ReadArraywInterp(solver->ndims,
99  adv_nvar,
100  solver->dim_global,
101  solver->dim_local,
102  dim_adv,
103  solver->ghosts,
104  solver,
105  mpi,
106  NULL,
107  adv,
108  fname_root,
109  &flag);
110  if (ierr) {
111  fprintf(stderr,"Error in LinearADRAdvectionField()\n");
112  fprintf(stderr," ReadArraywInterp() returned error!\n");
113  return ierr;
114  }
115  }
116 
117  if (!flag) {
118  /* if advection file not available, set it to zero */
119  _ArraySetValue_(adv,solver->npoints_local_wghosts*adv_nvar,0.0);
120  }
121 
122  /* if parallel, exchange MPI-boundary ghost point data */
124  adv_nvar,
125  solver->dim_local,
126  solver->ghosts,
127  mpi,
128  adv );
129 
130  /* Along each dimension:
131  - If boundaries are periodic, then set ghost points accordingly if
132  number of MPI ranks along that dimension is 1 (if more than 1,
133  MPIExchangeBoundariesnD() called above has taken care of it).
134  - Else, extrapolate at the boundaries. */
135  int indexb[solver->ndims],
136  indexi[solver->ndims],
137  bounds[solver->ndims],
138  offset[solver->ndims];
139  for (d = 0; d < solver->ndims; d++) {
140  if (solver->isPeriodic[d] && (mpi->iproc[d] == 1)) {
141  _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
142  /* left boundary */
143  done = 0; _ArraySetValue_(indexb,solver->ndims,0);
144  _ArraySetValue_(offset,solver->ndims,0); offset[d] = -ghosts;
145  while (!done) {
146  _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = indexb[d] + dim[d] - ghosts;
147  int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
148  int p2; _ArrayIndex1D_ (solver->ndims,dim,indexi,ghosts,p2);
149  _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
150  _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
151  }
152  /* right boundary */
153  done = 0; _ArraySetValue_(indexb,solver->ndims,0);
154  _ArraySetValue_(offset,solver->ndims,0); offset[d] = dim[d];
155  while (!done) {
156  _ArrayCopy1D_(indexb,indexi,solver->ndims);
157  int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
158  int p2; _ArrayIndex1D_ (solver->ndims,dim,indexi,ghosts,p2);
159  _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
160  _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
161  }
162  } else {
163  /* left boundary */
164  if (!mpi->ip[d]) {
165  _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
166  _ArraySetValue_(offset,solver->ndims,0); offset[d] = -ghosts;
167  done = 0; _ArraySetValue_(indexb,solver->ndims,0);
168  while (!done) {
169  _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = ghosts-1-indexb[d];
170  int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
171  int p2; _ArrayIndex1D_ (solver->ndims,dim,indexi,ghosts,p2);
172  _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
173  _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
174  }
175  }
176  /* right boundary */
177  if (mpi->ip[d] == mpi->iproc[d]-1) {
178  _ArrayCopy1D_(dim,bounds,solver->ndims); bounds[d] = ghosts;
179  _ArraySetValue_(offset,solver->ndims,0); offset[d] = dim[d];
180  done = 0; _ArraySetValue_(indexb,solver->ndims,0);
181  while (!done) {
182  _ArrayCopy1D_(indexb,indexi,solver->ndims); indexi[d] = dim[d]-1-indexb[d];
183  int p1; _ArrayIndex1DWO_(solver->ndims,dim,indexb,offset,ghosts,p1);
184  int p2; _ArrayIndex1D_ (solver->ndims,dim,indexi,ghosts,p2);
185  _ArrayCopy1D_((adv+p2*adv_nvar), (adv+p1*adv_nvar), adv_nvar);
186  _ArrayIncrementIndex_(solver->ndims,bounds,indexb,done);
187  }
188  }
189  }
190  }
191 
192  return(0);
193 }
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
char adv_filename[_MAX_STRING_SIZE_]
Definition: linearadr.h:43
int * dim_global
Definition: hypar.h:33
#define _ArrayIncrementIndex_(N, imax, i, done)
void GetStringFromInteger(int, char *, int)
int adv_arr_size
Definition: linearadr.h:47
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
double * a
Definition: linearadr.h:50
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
int ReadArray(int, int, int *, int *, int, void *, void *, double *, double *, char *, int *)
Definition: ReadArray.c:25
#define _MAX_STRING_SIZE_
Definition: basic.h:14
int * isPeriodic
Definition: hypar.h:162
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int ReadArraywInterp(int, int, int *, int *, int *, int, void *, void *, double *, double *, char *, int *)
int constant_advection
Definition: linearadr.h:40
int ndims
Definition: hypar.h:26
#define IERR
Definition: basic.h:16
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
double LinearADRComputeCFL ( void *  s,
void *  m,
double  dt,
double  t 
)

Computes the maximum CFL number over the domain. Note that the CFL is computed over the local domain on this processor only.

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables
dtTime step size for which to compute the CFL
tTime

Definition at line 15 of file LinearADRComputeCFL.c.

20 {
21  HyPar *solver = (HyPar*) s;
22  LinearADR *params = (LinearADR*) solver->physics;
23 
24  int ndims = solver->ndims;
25  int nvars = solver->nvars;
26  int ghosts = solver->ghosts;
27  int *dim = solver->dim_local;
28 
29  double max_cfl = 0;
30  if (params->constant_advection == 1) {
31 
32  int d, i, v;
33  for (d = 0; d < ndims; d++) {
34  for (i = 0; i < dim[d]; i++) {
35  for (v = 0; v < nvars; v++) {
36  double dxinv; _GetCoordinate_(d,i,dim,ghosts,solver->dxinv,dxinv);
37  double local_cfl = params->a[nvars*d+v]*dt*dxinv;
38  if (local_cfl > max_cfl) max_cfl = local_cfl;
39  }
40  }
41  }
42 
43  } else if (params->constant_advection == 0) {
44 
45  int d;
46  for (d = 0; d < ndims; d++) {
47  int index[ndims];
48  int done = 0; _ArraySetValue_(index,ndims,0);
49  while (!done) {
50  int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
51  double dxinv; _GetCoordinate_(0,index[0],dim,ghosts,solver->dxinv,dxinv);
52  int v;
53  for (v = 0; v < nvars; v++) {
54  double a = params->a[nvars*ndims*p+nvars*d+v];
55  double local_cfl = a*dt*dxinv;
56  if (local_cfl > max_cfl) max_cfl = local_cfl;
57  }
58  _ArrayIncrementIndex_(ndims,dim,index,done);
59  }
60  }
61 
62  }
63 
64  return(max_cfl);
65 }
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
double * a
Definition: linearadr.h:50
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
int nvars
Definition: hypar.h:29
double * dxinv
Definition: hypar.h:110
int constant_advection
Definition: linearadr.h:40
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double LinearADRComputeDiffNumber ( void *  s,
void *  m,
double  dt,
double  t 
)

Computes the maximum diffusion number over the domain. Note that the diffusion number is computed over the local domain on this processor only.

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables
dtTime step size for which to compute the CFL
tTime

Definition at line 14 of file LinearADRComputeDiffNumber.c.

19 {
20  HyPar *solver = (HyPar*) s;
21  LinearADR *params = (LinearADR*) solver->physics;
22  int d, i, v;
23 
24  int ndims = solver->ndims;
25  int nvars = solver->nvars;
26  int ghosts = solver->ghosts;
27  int *dim = solver->dim_local;
28 
29  double max_diffno = 0;
30  for (d = 0; d < ndims; d++) {
31  for (i = 0; i < dim[d]; i++) {
32  for (v = 0; v < nvars; v++) {
33  double dxinv; _GetCoordinate_(d,i,dim,ghosts,solver->dxinv,dxinv);
34  double local_diffno = params->d[nvars*d+v] * dt * dxinv * dxinv;
35  if (local_diffno > max_diffno) max_diffno = local_diffno;
36  }
37  }
38  }
39 
40  return(max_diffno);
41 }
double * d
Definition: linearadr.h:53
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
int ghosts
Definition: hypar.h:52
int nvars
Definition: hypar.h:29
double * dxinv
Definition: hypar.h:110
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int LinearADRAdvection ( double *  f,
double *  u,
int  dir,
void *  s,
double  t 
)

Evaluate the advection term in the linear advection-diffusion-reaction model:
Compute

\begin{equation} a_d u \end{equation}

given \(u\) and \(d\) in the hyperbolic term

\begin{equation} \sum_d \frac {\partial} {\partial x_d} \left( a_d u \right) \end{equation}

Parameters
fArray to hold the computed flux (same size and layout as u)
uArray containing the solution
dirSpatial dimension \(d\)
sSolver object of type HyPar
tCurrent time

Definition at line 37 of file LinearADRAdvection.c.

43 {
44  HyPar *solver = (HyPar*) s;
45  LinearADR *param = (LinearADR*) solver->physics;
46  double *adv = param->a;
47  int i, v;
48 
49  int *dim = solver->dim_local;
50  int ghosts = solver->ghosts;
51  int ndims = solver->ndims;
52  int nvars = solver->nvars;
53 
54  int index[ndims],
55  bounds[ndims],
56  offset[ndims];
57 
58  /* set bounds for array index to include ghost points */
59  _ArrayCopy1D_(dim,bounds,ndims);
60  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
61 
62  /* set offset such that index is compatible with ghost point arrangement */
63  _ArraySetValue_(offset,ndims,-ghosts);
64 
65  int done = 0; _ArraySetValue_(index,ndims,0);
66  if (param->constant_advection == 1) {
67  while (!done) {
68  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
69  for (v = 0; v < nvars; v++) {
70  f[nvars*p+v] = param->a[nvars*dir+v] * u[nvars*p+v];
71  }
72  _ArrayIncrementIndex_(ndims,bounds,index,done);
73  }
74  } else if (param->constant_advection == 0) {
75  while (!done) {
76  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
77  for (v = 0; v < nvars; v++) {
78  double a = adv[nvars*ndims*p+nvars*dir+v];
79  f[nvars*p+v] = a * u[nvars*p+v];
80  }
81  _ArrayIncrementIndex_(ndims,bounds,index,done);
82  }
83  } else {
84  while (!done) {
85  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
86  for (v = 0; v < nvars; v++) {
87  f[nvars*p+v] = 0.0;
88  }
89  _ArrayIncrementIndex_(ndims,bounds,index,done);
90  }
91  }
92 
93  return(0);
94 }
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
double * a
Definition: linearadr.h:50
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int constant_advection
Definition: linearadr.h:40
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int LinearADRDiffusionG ( double *  f,
double *  u,
int  dir,
void *  s,
double  t 
)

Evaluate the diffusion term in the linear advection-diffusion-reaction model for a "pure Laplacian" type operator (no cross derivatives):
Compute

\begin{equation} \nu_d u \end{equation}

given \(u\) and \(d\) in the parabolic term

\begin{equation} \sum_d \frac {\partial^2} {\partial x_d^2} \left( \nu_d u \right) \end{equation}

Parameters
fArray to hold the computed diffusion term (same size and layout as u)
uArray containing the solution
dirSpatial dimension (unused since this is a 1D system)
sSolver object of type HyPar
tCurrent time

Definition at line 26 of file LinearADRDiffusion.c.

32 {
33  HyPar *solver = (HyPar*) s;
34  LinearADR *param = (LinearADR*) solver->physics;
35  int i, v;
36 
37  int *dim = solver->dim_local;
38  int ghosts = solver->ghosts;
39  int ndims = solver->ndims;
40  int nvars = solver->nvars;
41  int index[ndims], bounds[ndims], offset[ndims];
42 
43  /* set bounds for array index to include ghost points */
44  _ArrayCopy1D_(dim,bounds,ndims);
45  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
46 
47  /* set offset such that index is compatible with ghost point arrangement */
48  _ArraySetValue_(offset,ndims,-ghosts);
49 
50  int done = 0; _ArraySetValue_(index,ndims,0);
51  while (!done) {
52  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
53  for (v = 0; v < nvars; v++) f[nvars*p+v] = param->d[nvars*dir+v] * u[nvars*p+v];
54  _ArrayIncrementIndex_(ndims,bounds,index,done);
55  }
56 
57  return(0);
58 }
#define _ArraySetValue_(x, size, value)
double * d
Definition: linearadr.h:53
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int LinearADRDiffusionH ( double *  f,
double *  u,
int  dir1,
int  dir2,
void *  s,
double  t 
)

Evaluate the diffusion term in the linear advection-diffusion-reaction model for a parabolic operator with no cross derivatives:
Compute

\begin{equation} \nu_d u \end{equation}

given \(u\) and \(d_1,d_2\) in the parabolic term

\begin{equation} \sum_{d_1}\sum_{d_2} \frac {\partial^2} {\partial x_{d_1} \partial x_{d_2}} \left( \nu_d u \right) \end{equation}

Note: it's not correctly implemented. Will implement when necessary.

Parameters
fArray to hold the computed diffusion term (same size and layout as u)
uArray containing the solution
dir1First spatial dimension of the double derivative \(d_1\)
dir2Second spatial dimension of the double derivative \(d_2\)
sSolver object of type HyPar
tCurrent time

Definition at line 74 of file LinearADRDiffusion.c.

81 {
82  HyPar *solver = (HyPar*) s;
83  LinearADR *param = (LinearADR*) solver->physics;
84  int i, v;
85 
86  int *dim = solver->dim_local;
87  int ghosts = solver->ghosts;
88  int ndims = solver->ndims;
89  int nvars = solver->nvars;
90  int index[ndims], bounds[ndims], offset[ndims];
91 
92  if (dir1 == dir2) {
93 
94  int dir = dir1;
95 
96  /* set bounds for array index to include ghost points */
97  _ArrayCopy1D_(dim,bounds,ndims);
98  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
99 
100  /* set offset such that index is compatible with ghost point arrangement */
101  _ArraySetValue_(offset,ndims,-ghosts);
102 
103  int done = 0; _ArraySetValue_(index,ndims,0);
104  while (!done) {
105  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
106  for (v = 0; v < nvars; v++) f[nvars*p+v] = param->d[nvars*dir+v] * u[nvars*p+v];
107  _ArrayIncrementIndex_(ndims,bounds,index,done);
108  }
109 
110  } else _ArraySetValue_(f,solver->npoints_local_wghosts*nvars,0.0);
111 
112 
113  return(0);
114 }
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
double * d
Definition: linearadr.h:53
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int LinearADRReaction ( )
int LinearADRUpwind ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Upwinding scheme for linear advection

Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension
sSolver object of type HyPar
tCurrent solution time

Definition at line 16 of file LinearADRUpwind.c.

26 {
27  HyPar *solver = (HyPar*) s;
28  LinearADR *param = (LinearADR*) solver->physics;
29  int done,v;
30 
31  int ndims = solver->ndims;
32  int nvars = solver->nvars;
33  int ghosts= solver->ghosts;
34  int *dim = solver->dim_local;
35 
36  double *a = param->a;
37 
38  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
39  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
40  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
41 
42  if (param->constant_advection == 1) {
43 
44  done = 0; _ArraySetValue_(index_outer,ndims,0);
45  while (!done) {
46  _ArrayCopy1D_(index_outer,index_inter,ndims);
47  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
48  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
49  for (v = 0; v < nvars; v++) {
50  fI[nvars*p+v] = (a[nvars*dir+v] > 0 ? fL[nvars*p+v] : fR[nvars*p+v] );
51  }
52  }
53  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
54  }
55 
56  } else if (param->constant_advection == 0) {
57 
58  done = 0; _ArraySetValue_(index_outer,ndims,0);
59  while (!done) {
60  _ArrayCopy1D_(index_outer,index_inter,ndims);
61  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
62  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
63  int indexL[ndims]; _ArrayCopy1D_(index_inter,indexL,ndims); indexL[dir]--;
64  int indexR[ndims]; _ArrayCopy1D_(index_inter,indexR,ndims);
65  int pL; _ArrayIndex1D_(ndims,dim,indexL,ghosts,pL);
66  int pR; _ArrayIndex1D_(ndims,dim,indexR,ghosts,pR);
67  for (v = 0; v < nvars; v++) {
68  double eigL = a[nvars*ndims*pL+nvars*dir+v],
69  eigR = a[nvars*ndims*pR+nvars*dir+v];
70  if ((eigL > 0) && (eigR > 0)) {
71  fI[nvars*p+v] = fL[nvars*p+v];
72  } else if ((eigL < 0) && (eigR < 0)) {
73  fI[nvars*p+v] = fR[nvars*p+v];
74  } else {
75  double alpha = max(absolute(eigL), absolute(eigR));
76  fI[nvars*p+v] = 0.5 * (fL[nvars*p+v] + fR[nvars*p+v] - alpha * (uR[nvars*p+v] - uL[nvars*p+v]));
77  }
78  }
79  }
80  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
81  }
82 
83  } else {
84 
85  done = 0; _ArraySetValue_(index_outer,ndims,0);
86  while (!done) {
87  _ArrayCopy1D_(index_outer,index_inter,ndims);
88  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
89  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
90  for (v = 0; v < nvars; v++) {
91  fI[nvars*p+v] = 0.0;
92  }
93  }
94  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
95  }
96 
97  }
98 
99  return(0);
100 }
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
double * a
Definition: linearadr.h:50
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int constant_advection
Definition: linearadr.h:40
int ndims
Definition: hypar.h:26
#define absolute(a)
Definition: math_ops.h:32
#define max(a, b)
Definition: math_ops.h:18
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int LinearADRCenteredFlux ( double *  fI,
double *  fL,
double *  fR,
double *  uL,
double *  uR,
double *  u,
int  dir,
void *  s,
double  t 
)

Centered scheme for linear advection

Parameters
fIComputed upwind interface flux
fLLeft-biased reconstructed interface flux
fRRight-biased reconstructed interface flux
uLLeft-biased reconstructed interface solution
uRRight-biased reconstructed interface solution
uCell-centered solution
dirSpatial dimension
sSolver object of type HyPar
tCurrent solution time

Definition at line 103 of file LinearADRUpwind.c.

113 {
114  HyPar *solver = (HyPar*) s;
115  LinearADR *param = (LinearADR*) solver->physics;
116 
117  int ndims = solver->ndims;
118  int nvars = solver->nvars;
119  int ghosts= solver->ghosts;
120  int *dim = solver->dim_local;
121 
122  int index_outer[ndims], index_inter[ndims], bounds_outer[ndims], bounds_inter[ndims];
123  _ArrayCopy1D_(dim,bounds_outer,ndims); bounds_outer[dir] = 1;
124  _ArrayCopy1D_(dim,bounds_inter,ndims); bounds_inter[dir] += 1;
125 
126  int done = 0; _ArraySetValue_(index_outer,ndims,0);
127  while (!done) {
128  _ArrayCopy1D_(index_outer,index_inter,ndims);
129  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
130  int p; _ArrayIndex1D_(ndims,bounds_inter,index_inter,0,p);
131  for (int v = 0; v < nvars; v++) {
132  fI[nvars*p+v] = 0.5 * (fL[nvars*p+v] + fR[nvars*p+v]);
133  }
134  }
135  _ArrayIncrementIndex_(ndims,bounds_outer,index_outer,done);
136  }
137 
138  return(0);
139 }
#define _ArraySetValue_(x, size, value)
#define _ArrayIncrementIndex_(N, imax, i, done)
void * physics
Definition: hypar.h:266
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayCopy1D_(x, y, size)
int nvars
Definition: hypar.h:29
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int LinearADRWriteAdvField ( void *  s,
void *  m,
double  a_t 
)

Write out the advection field to file

Parameters
sSolver object of type HyPar
mMPI object of type MPIVariables
a_tCurrent simulation time

Definition at line 18 of file LinearADRWriteAdvField.c.

21 {
22  HyPar *solver = (HyPar*) s;
23  MPIVariables *mpi = (MPIVariables*) m;
24  LinearADR *param = (LinearADR*) solver->physics;
25 
26  if (param->constant_advection == 0) {
27 
28  char fname_root[_MAX_STRING_SIZE_] = "advection_field";
29  if (solver->nsims > 1) {
30  char index[_MAX_STRING_SIZE_];
31  GetStringFromInteger(solver->my_idx, index, (int)log10(solver->nsims)+1);
32  strcat(fname_root, "_");
33  strcat(fname_root, index);
34  }
35 
36  int adv_nvar = solver->ndims * solver->nvars;
37  WriteArray( solver->ndims,
38  adv_nvar,
39  solver->dim_global,
40  solver->dim_local,
41  solver->ghosts,
42  solver->x,
43  param->a,
44  solver,mpi,
45  fname_root );
46 
47  if (!strcmp(solver->plot_solution, "yes")) {
48  PlotArray( solver->ndims,
49  adv_nvar,
50  solver->dim_global,
51  solver->dim_local,
52  solver->ghosts,
53  solver->x,
54  param->a,
55  a_t,
56  solver,mpi,
57  fname_root );
58  }
59 
60  }
61 
62  return(0);
63 }
int * dim_global
Definition: hypar.h:33
void GetStringFromInteger(int, char *, int)
void * physics
Definition: hypar.h:266
char plot_solution[_MAX_STRING_SIZE_]
Definition: hypar.h:194
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int * dim_local
Definition: hypar.h:37
double * a
Definition: linearadr.h:50
int ghosts
Definition: hypar.h:52
int nsims
Definition: hypar.h:64
int WriteArray(int, int, int *, int *, int, double *, double *, void *, void *, char *)
Definition: WriteArray.c:27
#define _MAX_STRING_SIZE_
Definition: basic.h:14
int nvars
Definition: hypar.h:29
int PlotArray(int, int, int *, int *, int, double *, double *, double, void *, void *, char *)
Definition: PlotArray.c:31
int constant_advection
Definition: linearadr.h:40
int ndims
Definition: hypar.h:26
int my_idx
Definition: hypar.h:61
Structure of MPI-related variables.
double * x
Definition: hypar.h:107
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int LinearADRAdvectionJacobian ( double *  Jac,
double *  u,
void *  p,
int  dir,
int  nvars,
int  upw 
)

Function to compute the flux Jacobian for the hyperbolic (advection) part of the linear-advection-diffusion-reaction model.

Parameters
JacJacobian matrix of size 1 (nvar = 1)
usolution at a grid point
pobject containing physics-related parameters
dirdimension (x/y/z)
nvarsnumber of components
upw0 -> send back complete Jacobian, 1 -> send back Jacobian of right(+)-moving flux, -1 -> send back Jacobian of left(-)-moving flux

Definition at line 11 of file LinearADRJacobian.c.

19 {
20  LinearADR *param = (LinearADR*) p;
21 
22  if (param->a) {
23  *Jac = (1-absolute(upw))*absolute(param->a[dir])
24  + absolute(upw) * (1+upw) * max(0,param->a[dir]) * 0.5
25  - absolute(upw) * (1-upw) * min(0,param->a[dir]) * 0.5 ;
26  } else {
27  /* no advection term */
28  *Jac = 0.0;
29  }
30 
31  return 0;
32 }
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
double * a
Definition: linearadr.h:50
#define absolute(a)
Definition: math_ops.h:32
#define max(a, b)
Definition: math_ops.h:18
#define min(a, b)
Definition: math_ops.h:14
int LinearADRDiffusionJacobian ( double *  Jac,
double *  u,
void *  p,
int  dir,
int  nvars 
)

Function to compute the Jacobian for the parabolic (diffusion) part of the linear-advection-diffusion-reaction model.

Parameters
JacJacobian matrix of size 1 (nvar = 1)
usolution at a grid point
pobject containing physics-related parameters
dirdimension (x/y/z)
nvarsnumber of components

Definition at line 36 of file LinearADRJacobian.c.

41 {
42  LinearADR *param = (LinearADR*) p;
43 
44  int v;
45  for (v = 0; v < nvars; v++) {
46  Jac[nvars*v+v] = -param->d[nvars*dir+v];
47  }
48 
49  return 0;
50 }
double * d
Definition: linearadr.h:53
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int LinearADRInitialize ( void *  s,
void *  m 
)

Initialize the linear advection-diffusion-reaction physics module - allocate and set physics-related parameters, read physics-related inputs from file, and set the physics-related function pointers in HyPar

This file reads the file "physics.inp" that must have the following format:

begin
    <keyword>   <value>
    <keyword>   <value>
    <keyword>   <value>
    ...
    <keyword>   <value>
end

where the list of keywords are:

Keyword name Type Variable Default value
advection_filename char[] LinearADR::adv_filename "none"
advection double[] LinearADR::a 0
diffusion double[] LinearADR::d 0
centered_flux char[] LinearADR::centered_flux "no"

Note:

  • "physics.inp" is optional; if absent, default values will be used.
  • Please do not specify both "advection" and "advection_filename"!
Parameters
sSolver object of type HyPar
mObject of type MPIVariables containing MPI-related info

Definition at line 59 of file LinearADRInitialize.c.

61 {
62  HyPar *solver = (HyPar*) s;
63  MPIVariables *mpi = (MPIVariables*) m;
64  LinearADR *physics = (LinearADR*) solver->physics;
65  int i,ferr;
66 
67  static int count = 0;
68 
69  /* default values */
70  physics->constant_advection = -1;
71  strcpy(physics->adv_filename,"none");
72  physics->a = NULL;
73  physics->adv_arr_size = -1;
74  strcpy(physics->centered_flux,"no");
75 
76  physics->d = (double*) calloc (solver->ndims*solver->nvars,sizeof(double));
77  _ArraySetValue_(physics->d,solver->ndims*solver->nvars,0.0);
78 
79  /* reading physical model specific inputs - all processes */
80  if (!mpi->rank) {
81 
82  FILE *in;
83  if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
84  in = fopen("physics.inp","r");
85 
86  if (!in) {
87 
88  fprintf(stderr,"Error: File \"physics.inp\" not found.\n");
89  return(1);
90 
91  } else {
92 
93  char word[_MAX_STRING_SIZE_];
94  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
95  if (!strcmp(word, "begin")) {
96  while (strcmp(word, "end")) {
97 
98  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
99 
100  if (!strcmp(word, "advection_filename")) {
101  if (physics->constant_advection != -1) {
102  fprintf(stderr,"Error in LinearADRInitialize():\n");
103  fprintf(stderr,"Maybe advection_filename and advection are both specified.\n");
104  return 1;
105  } else {
106  /* read filename for spatially-varying advection field */
107  physics->constant_advection = 0;
108  physics->adv_arr_size = solver->npoints_local_wghosts*solver->ndims*solver->nvars;
109  physics->a = (double*) calloc ( physics->adv_arr_size, sizeof(double) );
110  _ArraySetValue_(physics->a, physics->adv_arr_size, 0.0);
111  ferr = fscanf(in,"%s",physics->adv_filename); if (ferr != 1) return(ferr);
112  }
113  } else if (!strcmp(word, "advection")) {
114  if (physics->constant_advection != -1) {
115  fprintf(stderr,"Error in LinearADRInitialize():\n");
116  fprintf(stderr,"Maybe advection_filename and advection are both specified.\n");
117  return 1;
118  } else {
119  /* read advection coefficients */
120  physics->constant_advection = 1;
121  physics->adv_arr_size = solver->ndims*solver->nvars;
122  physics->a = (double*) calloc ( physics->adv_arr_size, sizeof(double) );
123  for (i=0; i<solver->ndims*solver->nvars; i++) ferr = fscanf(in,"%lf",&physics->a[i]);
124  if (ferr != 1) return(1);
125  }
126  } else if (!strcmp(word, "diffusion")) {
127  /* read diffusion coefficients */
128  for (i=0; i<solver->ndims*solver->nvars; i++) ferr = fscanf(in,"%lf",&physics->d[i]);
129  if (ferr != 1) return(1);
130  } else if (!strcmp(word, "centered_flux")) {
131  ferr = fscanf(in, "%s", physics->centered_flux);
132  if (ferr != 1) return(ferr);
133  } else if (strcmp(word,"end")) {
134  char useless[_MAX_STRING_SIZE_];
135  ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
136  printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",word,useless);
137  printf("recognized or extraneous. Ignoring.\n");
138  }
139  }
140 
141  } else {
142 
143  fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
144  return(1);
145 
146  }
147  }
148 
149  fclose(in);
150 
151  }
152 
153 #ifndef serial
154  MPIBroadcast_integer(&physics->constant_advection,1,0,&mpi->world);
155  MPIBroadcast_integer(&physics->adv_arr_size,1,0,&mpi->world);
156 #endif
157 
158  if (mpi->rank) {
159  physics->a = (double*) calloc (physics->adv_arr_size,sizeof(double));
160  _ArraySetValue_(physics->a, physics->adv_arr_size, 0.0);
161  }
162 
163  if (physics->constant_advection == 1) {
164 #ifndef serial
165  MPIBroadcast_double(physics->a,physics->adv_arr_size,0,&mpi->world);
166 #endif
167  } else if (physics->constant_advection == 0) {
168 #ifndef serial
170 #endif
171  }
172 
173 #ifndef serial
174  MPIBroadcast_double(physics->d,solver->ndims*solver->nvars,0,&mpi->world);
176 #endif
177 
178  if (!strcmp(solver->SplitHyperbolicFlux,"yes")) {
179  if (!mpi->rank) {
180  fprintf(stderr,"Error in LinearADRInitialize: This physical model does not have a splitting ");
181  fprintf(stderr,"of the hyperbolic term defined.\n");
182  }
183  return(1);
184  }
185 
186  /* initializing physical model-specific functions */
189  solver->FFunction = LinearADRAdvection;
190  solver->GFunction = LinearADRDiffusionG;
191  solver->HFunction = LinearADRDiffusionH;
192  solver->SFunction = LinearADRReaction;
195 
196  if (!strcmp(physics->centered_flux,"no")) {
197  solver->Upwind = LinearADRUpwind;
198  } else {
199  solver->Upwind = LinearADRCenteredFlux;
200  }
201 
202  if (physics->constant_advection == 0) {
205  } else {
206  solver->PhysicsInput = NULL;
207  solver->PhysicsOutput = NULL;
208  }
209 
210  count++;
211  return(0);
212 }
int npoints_local_wghosts
Definition: hypar.h:42
#define _ArraySetValue_(x, size, value)
int(* JFunction)(double *, double *, void *, int, int, int)
Definition: hypar.h:326
double(* ComputeCFL)(void *, void *, double, double)
Definition: hypar.h:269
char adv_filename[_MAX_STRING_SIZE_]
Definition: linearadr.h:43
double * d
Definition: linearadr.h:53
int MPIBroadcast_double(double *, int, int, void *)
Definition: MPIBroadcast.c:9
int adv_arr_size
Definition: linearadr.h:47
int(* HFunction)(double *, double *, int, int, void *, double)
Definition: hypar.h:313
void * physics
Definition: hypar.h:266
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
Structure containing variables and parameters specific to the linear advection-diffusion-reaction mod...
Definition: linearadr.h:37
int(* PhysicsOutput)(void *, void *, double)
Definition: hypar.h:347
double * a
Definition: linearadr.h:50
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int(* PhysicsInput)(void *, void *, int, int, int *)
Definition: hypar.h:351
int LinearADRWriteAdvField(void *, void *, double)
int LinearADRCenteredFlux(double *, double *, double *, double *, double *, double *, int, void *, double)
MPI_Comm world
#define _MAX_STRING_SIZE_
Definition: basic.h:14
double LinearADRComputeCFL(void *s, void *m, double dt, double t)
int LinearADRDiffusionG(double *f, double *u, int dir, void *s, double t)
double(* ComputeDiffNumber)(void *, void *, double, double)
Definition: hypar.h:272
int LinearADRDiffusionJacobian(double *, double *, void *, int, int)
int LinearADRAdvection(double *f, double *u, int dir, void *s, double t)
int LinearADRDiffusionH(double *f, double *u, int dir1, int dir2, void *s, double t)
int LinearADRAdvectionField(void *s, void *m, int idx, int nsims, int *dim_adv)
int(* SFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:317
int(* KFunction)(double *, double *, void *, int, int)
Definition: hypar.h:331
int nvars
Definition: hypar.h:29
char centered_flux[_MAX_STRING_SIZE_]
Definition: linearadr.h:57
int MPIBroadcast_character(char *, int, int, void *)
Definition: MPIBroadcast.c:37
int constant_advection
Definition: linearadr.h:40
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
int ndims
Definition: hypar.h:26
int LinearADRReaction()
int LinearADRAdvectionJacobian(double *, double *, void *, int, int, int)
double LinearADRComputeDiffNumber(void *s, void *m, double dt, double t)
Structure of MPI-related variables.
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
int LinearADRUpwind(double *, double *, double *, double *, double *, double *, int, void *, double)
int(* GFunction)(double *, double *, int, void *, double)
Definition: hypar.h:310
int MPIBroadcast_integer(int *, int, int, void *)
Definition: MPIBroadcast.c:23