HyPar  1.0
Finite-Difference Hyperbolic-Parabolic PDE Solver on Cartesian Grids
Numa2DInitialize.c File Reference
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <basic.h>
#include <arrayfunctions.h>
#include <mathfunctions.h>
#include <boundaryconditions.h>
#include <physicalmodels/numa2d.h>
#include <mpivars.h>
#include <hypar.h>

Go to the source code of this file.

Functions

double Numa2DComputeCFL (void *, void *, double, double)
 
int Numa2DFlux (double *, double *, int, void *, double)
 
int Numa2DStiffFlux (double *, double *, int, void *, double)
 
int Numa2DSource (double *, double *, void *, void *, double)
 
int Numa2DParabolicFunction (double *, double *, void *, void *, double)
 
int Numa2DRusanovFlux (double *, double *, double *, double *, double *, double *, int, void *, double)
 
int Numa2DRusanovLinearFlux (double *, double *, double *, double *, double *, double *, int, void *, double)
 
void Numa2DCalculateStandardAtmosphere_1 (void *, double, double *, double *, double *, double *)
 
void Numa2DCalculateStandardAtmosphere_2 (void *, double, double *, double *, double *, double *)
 
int Numa2DInitialize (void *s, void *m)
 

Function Documentation

◆ Numa2DComputeCFL()

double Numa2DComputeCFL ( void *  ,
void *  ,
double  ,
double   
)

Definition at line 9 of file Numa2DComputeCFL.c.

10 {
11  HyPar *solver = (HyPar*) s;
12  Numa2D *param = (Numa2D*) solver->physics;
13 
14  int *dim = solver->dim_local;
15  int ghosts = solver->ghosts;
16  int ndims = solver->ndims;
17  double *u = solver->u;
18  int index[ndims];
19 
20  double max_cfl = 0;
21  int done = 0; _ArraySetValue_(index,ndims,0);
22  while (!done) {
23  int p; _ArrayIndex1D_(ndims,dim,index,ghosts,p);
24  double drho,uvel,vvel,dT,rho0,T0,P0,EP,c,ycoord;
25 
26  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->x,ycoord);
27  param->StandardAtmosphere (param,ycoord,&EP,&P0,&rho0,&T0);
28  _Numa2DGetFlowVars_ ((u+_MODEL_NVARS_*p),drho,uvel,vvel,dT,rho0);
29  _Numa2DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
30 
31  double dxinv, dyinv;
32  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv); /* 1/dx */
33  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv); /* 1/dy */
34 
35  double local_cfl[2];
36  local_cfl[_XDIR_] = (absolute(uvel)+c)*dt*dxinv; /* local cfl for this grid point (x) */
37  local_cfl[_YDIR_] = (absolute(vvel)+c)*dt*dyinv; /* local cfl for this grid point (y) */
38  if (local_cfl[_XDIR_] > max_cfl) max_cfl = local_cfl[_XDIR_];
39  if (local_cfl[_YDIR_] > max_cfl) max_cfl = local_cfl[_YDIR_];
40 
41  _ArrayIncrementIndex_(ndims,dim,index,done);
42  }
43 
44  return(max_cfl);
45 }
#define _Numa2DComputeSpeedofSound_(gamma, R, T0, dT, rho0, drho, EP, c)
Definition: numa2d.h:99
#define absolute(a)
Definition: math_ops.h:32
double * u
Definition: hypar.h:116
double * x
Definition: hypar.h:107
double R
Definition: numa2d.h:111
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Definition: numa2d.h:109
#define _YDIR_
Definition: euler2d.h:41
double gamma
Definition: numa2d.h:110
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _Numa2DGetFlowVars_(u, drho, uvel, vvel, dT, rho0)
Definition: numa2d.h:33
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
double * dxinv
Definition: hypar.h:110

◆ Numa2DFlux()

int Numa2DFlux ( double *  ,
double *  ,
int  ,
void *  ,
double   
)

Definition at line 7 of file Numa2DFlux.c.

8 {
9  HyPar *solver = (HyPar*) s;
10  Numa2D *param = (Numa2D*) solver->physics;
11  int i;
12 
13  int *dim = solver->dim_local;
14  int ghosts = solver->ghosts;
15  int ndims = solver->ndims;
16  int index[ndims], bounds[ndims], offset[ndims];
17 
18  /* set bounds for array index to include ghost points */
19  _ArrayCopy1D_(dim,bounds,ndims);
20  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
21 
22  /* set offset such that index is compatible with ghost point arrangement */
23  _ArraySetValue_(offset,ndims,-ghosts);
24 
25  int done = 0; _ArraySetValue_(index,ndims,0);
26  while (!done) {
27  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
28  double drho,uvel,vvel,dT,dP,rho0,T0,P0,EP,ycoord;
29 
30  _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
31  param->StandardAtmosphere(param,ycoord,&EP,&P0,&rho0,&T0);
32 
33  _Numa2DGetFlowVars_ ((u+_MODEL_NVARS_*p),drho,uvel,vvel,dT,rho0);
34  _Numa2DComputePressure_ (param,T0,dT,P0,dP);
35  _Numa2DSetFlux_ ((f+_MODEL_NVARS_*p),dir,drho,uvel,vvel,dT,dP,rho0,T0);
36 
37  _ArrayIncrementIndex_(ndims,bounds,index,done);
38  }
39 
40  return(0);
41 }
double * x
Definition: hypar.h:107
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Definition: numa2d.h:109
#define _YDIR_
Definition: euler2d.h:41
#define _Numa2DComputePressure_(params, T0, dT, P0, dP)
Definition: numa2d.h:79
#define _Numa2DSetFlux_(f, dir, drho, uvel, vvel, dT, dP, rho0, T0)
Definition: numa2d.h:41
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _Numa2DGetFlowVars_(u, drho, uvel, vvel, dT, rho0)
Definition: numa2d.h:33
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
#define _ArrayCopy1D_(x, y, size)

◆ Numa2DStiffFlux()

int Numa2DStiffFlux ( double *  ,
double *  ,
int  ,
void *  ,
double   
)

Definition at line 43 of file Numa2DFlux.c.

44 {
45  HyPar *solver = (HyPar*) s;
46  Numa2D *param = (Numa2D*) solver->physics;
47  int i;
48 
49  int *dim = solver->dim_local;
50  int ghosts = solver->ghosts;
51  int ndims = solver->ndims;
52  int index[ndims], bounds[ndims], offset[ndims];
53 
54  /* set bounds for array index to include ghost points */
55  _ArrayCopy1D_(dim,bounds,ndims);
56  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
57 
58  /* set offset such that index is compatible with ghost point arrangement */
59  _ArraySetValue_(offset,ndims,-ghosts);
60 
61  int done = 0; _ArraySetValue_(index,ndims,0);
62  while (!done) {
63  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
64  double drho,uvel,vvel,dT,dP,rho0,T0,P0,EP,ycoord;
65 
66  _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
67  param->StandardAtmosphere(param,ycoord,&EP,&P0,&rho0,&T0);
68 
69  _Numa2DGetFlowVars_ ((u+_MODEL_NVARS_*p),drho,uvel,vvel,dT,rho0);
70  _Numa2DComputeLinearizedPressure_ (param,T0,dT,P0,dP);
71  _Numa2DSetLinearFlux_ ((f+_MODEL_NVARS_*p),dir,drho,uvel,vvel,dT,dP,rho0,T0);
72 
73  _ArrayIncrementIndex_(ndims,bounds,index,done);
74  }
75 
76  return(0);
77 }
#define _Numa2DComputeLinearizedPressure_(params, T0, dT, P0, dP)
Definition: numa2d.h:88
double * x
Definition: hypar.h:107
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Definition: numa2d.h:109
#define _YDIR_
Definition: euler2d.h:41
#define _Numa2DSetLinearFlux_(f, dir, drho, uvel, vvel, dT, dP, rho0, T0)
Definition: numa2d.h:56
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _Numa2DGetFlowVars_(u, drho, uvel, vvel, dT, rho0)
Definition: numa2d.h:33
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
#define _ArrayCopy1D_(x, y, size)

◆ Numa2DSource()

int Numa2DSource ( double *  ,
double *  ,
void *  ,
void *  ,
double   
)

Definition at line 7 of file Numa2DSource.c.

8 {
9  HyPar *solver = (HyPar*) s;
10  Numa2D *param = (Numa2D*) solver->physics;
11  int i;
12 
13  int *dim = solver->dim_local;
14  int ghosts = solver->ghosts;
15  int ndims = solver->ndims;
16  int index[ndims], bounds[ndims], offset[ndims];
17 
18  /* set bounds for array index to include ghost points */
19  _ArrayCopy1D_(dim,bounds,ndims);
20  for (i=0; i<ndims; i++) bounds[i] += 2*ghosts;
21 
22  /* set offset such that index is compatible with ghost point arrangement */
23  _ArraySetValue_(offset,ndims,-ghosts);
24 
25  int done = 0; _ArraySetValue_(index,ndims,0);
26  while (!done) {
27  int p; _ArrayIndex1DWO_(ndims,dim,index,offset,ghosts,p);
28  double drho,uvel,vvel,dT,rho0,P0,EP,T0,ycoord;
29 
30  _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
31  param->StandardAtmosphere (param,ycoord,&EP,&P0,&rho0,&T0);
32  _Numa2DGetFlowVars_ ((u+_MODEL_NVARS_*p),drho,uvel,vvel,dT,rho0);
33  _Numa2DSetSource_ ((S+_MODEL_NVARS_*p),param,drho);
34 
35  /* some useless statements to avoid compiler warnings */
36  uvel = dT;
37  vvel = uvel;
38  dT = vvel;
39 
40  _ArrayIncrementIndex_(ndims,bounds,index,done);
41  }
42 
43  return(0);
44 }
#define _Numa2DSetSource_(s, param, drho)
Definition: numa2d.h:71
double * x
Definition: hypar.h:107
int ndims
Definition: hypar.h:26
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Definition: numa2d.h:109
#define _YDIR_
Definition: euler2d.h:41
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _Numa2DGetFlowVars_(u, drho, uvel, vvel, dT, rho0)
Definition: numa2d.h:33
int ghosts
Definition: hypar.h:52
#define _MODEL_NVARS_
Definition: euler1d.h:58
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
#define _ArrayCopy1D_(x, y, size)

◆ Numa2DParabolicFunction()

int Numa2DParabolicFunction ( double *  ,
double *  ,
void *  ,
void *  ,
double   
)

Definition at line 17 of file Numa2DParabolicFunction.c.

18 {
19  HyPar *solver = (HyPar*) s;
20  MPIVariables *mpi = (MPIVariables*) m;
21  Numa2D *physics = (Numa2D*) solver->physics;
22  int i,v,done;
23  double dxinv, dyinv;
25 
26  int ghosts = solver->ghosts;
27  int imax = solver->dim_local[0];
28  int jmax = solver->dim_local[1];
29  int *dim = solver->dim_local;
30  int size = (imax+2*ghosts)*(jmax+2*ghosts)*_MODEL_NVARS_;
31 
32  _ArraySetValue_(par,size,0.0);
33  if (physics->mu <= 0) return(0); /* inviscid flow */
34  solver->count_par++;
35 
36  double mu = physics->mu;
37 
38  /* allocate some arrays */
39  double *Q, *QDeriv, *FViscous, *FDeriv;
40  Q = (double*) calloc (size,sizeof(double)); /* primitive variables */
41  QDeriv = (double*) calloc (size,sizeof(double)); /* derivative of primitive variables */
42  FViscous = (double*) calloc (size,sizeof(double)); /* viscous flux */
43  FDeriv = (double*) calloc (size,sizeof(double)); /* derivative of viscous flux */
44 
45  int index[_MODEL_NDIMS_], bounds[_MODEL_NDIMS_], offset[_MODEL_NDIMS_];
46  /* set bounds for array index to include ghost points */
47  _ArrayCopy1D_(dim,bounds,_MODEL_NDIMS_); for (i=0; i<_MODEL_NDIMS_; i++) bounds[i] += 2*ghosts;
48  /* set offset such that index is compatible with ghost point arrangement */
49  _ArraySetValue_(offset,_MODEL_NDIMS_,-ghosts);
50 
51  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
52  while (!done) {
53  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
54  double drho,uvel,vvel,dT,rho0,T0,P0,EP,ycoord;
55 
56  _GetCoordinate_ (_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->x,ycoord);
57  physics->StandardAtmosphere (physics,ycoord,&EP,&P0,&rho0,&T0);
58  _Numa2DGetFlowVars_ ( (u+_MODEL_NVARS_*p),drho,uvel,vvel,dT,rho0);
59 
60  Q[_MODEL_NVARS_*p+0] = rho0 + drho; /* density */
61  Q[_MODEL_NVARS_*p+1] = uvel; /* x-velocity */
62  Q[_MODEL_NVARS_*p+2] = vvel; /* y-velocity */
63  Q[_MODEL_NVARS_*p+3] = (dT+T0)/(drho+rho0); /* potential temperature*/
64 
65  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
66  }
67 
68  /* Along X */
69  _ArraySetValue_(QDeriv,size,0.0);
70  IERR solver->FirstDerivativePar(QDeriv,Q,_XDIR_,1,solver,mpi); CHECKERR(ierr);
72  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
73  while (!done) {
74  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
75  _GetCoordinate_(_XDIR_,index[_XDIR_]-ghosts,dim,ghosts,solver->dxinv,dxinv);
76 
77  double rho, ux, vx, tx;
78  rho = Q [_MODEL_NVARS_*p+0];
79  ux = QDeriv[_MODEL_NVARS_*p+1] / dxinv;
80  vx = QDeriv[_MODEL_NVARS_*p+2] / dxinv;
81  tx = QDeriv[_MODEL_NVARS_*p+3] / dxinv;
82 
83  FViscous[_MODEL_NVARS_*p+0] = 0.0;
84  FViscous[_MODEL_NVARS_*p+1] = mu * rho * ux;
85  FViscous[_MODEL_NVARS_*p+2] = mu * rho * vx;
86  FViscous[_MODEL_NVARS_*p+3] = mu * rho * tx;
87 
88  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
89  }
90  _ArraySetValue_(FDeriv,size,0.0);
91  IERR solver->FirstDerivativePar(FDeriv,FViscous,_XDIR_,-1,solver,mpi); CHECKERR(ierr);
92  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
93  while (!done) {
94  int p; _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
95  _GetCoordinate_(_XDIR_,index[_XDIR_],dim,ghosts,solver->dxinv,dxinv);
96  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dxinv * (FDeriv+p)[v] );
97  _ArrayIncrementIndex_(_MODEL_NDIMS_,dim,index,done);
98  }
99 
100  /* Along Y */
101  _ArraySetValue_(QDeriv,size,0.0);
102  IERR solver->FirstDerivativePar(QDeriv,Q,_YDIR_,1,solver,mpi); CHECKERR(ierr);
103  IERR MPIExchangeBoundariesnD(_MODEL_NDIMS_,_MODEL_NVARS_,dim,ghosts,mpi,QDeriv);CHECKERR(ierr);
104  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
105  while (!done) {
106  int p; _ArrayIndex1DWO_(_MODEL_NDIMS_,dim,index,offset,ghosts,p);
107  _GetCoordinate_(_YDIR_,index[_YDIR_]-ghosts,dim,ghosts,solver->dxinv,dyinv);
108 
109  double rho, uy, vy, ty;
110  rho = Q [_MODEL_NVARS_*p+0];
111  uy = QDeriv[_MODEL_NVARS_*p+1] / dyinv;
112  vy = QDeriv[_MODEL_NVARS_*p+2] / dyinv;
113  ty = QDeriv[_MODEL_NVARS_*p+3] / dyinv;
114 
115  FViscous[_MODEL_NVARS_*p+0] = 0.0;
116  FViscous[_MODEL_NVARS_*p+1] = mu * rho * uy;
117  FViscous[_MODEL_NVARS_*p+2] = mu * rho * vy;
118  FViscous[_MODEL_NVARS_*p+3] = mu * rho * ty;
119 
120  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds,index,done);
121  }
122  _ArraySetValue_(FDeriv,size,0.0);
123  IERR solver->FirstDerivativePar(FDeriv,FViscous,_YDIR_,-1,solver,mpi); CHECKERR(ierr);
124  done = 0; _ArraySetValue_(index,_MODEL_NDIMS_,0);
125  while (!done) {
126  int p; _ArrayIndex1D_(_MODEL_NDIMS_,dim,index,ghosts,p); p *= _MODEL_NVARS_;
127  _GetCoordinate_(_YDIR_,index[_YDIR_],dim,ghosts,solver->dxinv,dyinv);
128  for (v=0; v<_MODEL_NVARS_; v++) (par+p)[v] += (dyinv * (FDeriv+p)[v] );
129  _ArrayIncrementIndex_(_MODEL_NDIMS_,dim,index,done);
130  }
131 
132  free(Q);
133  free(QDeriv);
134  free(FViscous);
135  free(FDeriv);
136 
137  return(0);
138 }
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
double mu
Definition: numa2d.h:114
int MPIExchangeBoundariesnD(int, int, int *, int, void *, double *)
double * x
Definition: hypar.h:107
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Definition: numa2d.h:109
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
int count_par
Definition: hypar.h:418
#define _ArrayIndex1D_(N, imax, i, ghost, index)
#define _ArrayIndex1DWO_(N, imax, i, offset, ghost, index)
#define _ArraySetValue_(x, size, value)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
int(* FirstDerivativePar)(double *, double *, int, int, void *, void *)
Definition: hypar.h:243
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _Numa2DGetFlowVars_(u, drho, uvel, vvel, dT, rho0)
Definition: numa2d.h:33
#define _XDIR_
Definition: euler1d.h:75
int ghosts
Definition: hypar.h:52
Structure of MPI-related variables.
#define _MODEL_NVARS_
Definition: euler1d.h:58
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
#define _DECLARE_IERR_
Definition: basic.h:17
#define _ArrayCopy1D_(x, y, size)
double * dxinv
Definition: hypar.h:110

◆ Numa2DRusanovFlux()

int Numa2DRusanovFlux ( double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
int  ,
void *  ,
double   
)

Definition at line 9 of file Numa2DUpwind.c.

10 {
11  HyPar *solver = (HyPar*) s;
12  Numa2D *param = (Numa2D*) solver->physics;
13  int done;
14 
15  int *dim = solver->dim_local;
16  int ghosts = solver->ghosts;
17 
18  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
19  _ArrayCopy1D2_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
20  _ArrayCopy1D2_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
21 
22  done = 0; int index_outer[2] = {0,0}, index_inter[2];
23  while (!done) {
24  _ArrayCopy1D2_(index_outer,index_inter,_MODEL_NDIMS_);
25  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
26  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
27  double udiff[_MODEL_NVARS_],drho,vel[2],dT,rho0,P0,T0,EP,c;
28  double ycoordL, ycoordR;
29 
30  if (dir == _YDIR_) {
31  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_]-1),dim,ghosts,solver->x,ycoordL);
32  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_] ),dim,ghosts,solver->x,ycoordR);
33  } else {
34  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_] ),dim,ghosts,solver->x,ycoordL);
35  ycoordR = ycoordL;
36  }
37 
38  /* Rusanov's upwinding scheme */
39 
40  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
41  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
42  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
43  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
44 
45  /* left of the interface */
46  param->StandardAtmosphere (param,ycoordL,&EP,&P0,&rho0,&T0);
47  _Numa2DGetFlowVars_ ((uL+_MODEL_NVARS_*p),drho,vel[0],vel[1],dT,rho0);
48  _Numa2DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
49  double alphaL = c + absolute(vel[dir]);
50 
51  /* right of the interface */
52  param->StandardAtmosphere (param,ycoordR,&EP,&P0,&rho0,&T0);
53  _Numa2DGetFlowVars_ ((uR+_MODEL_NVARS_*p),drho,vel[0],vel[1],dT,rho0);
54  _Numa2DComputeSpeedofSound_ (param->gamma,param->R,T0,dT,rho0,drho,EP,c);
55  double alphaR = c + absolute(vel[dir]);
56 
57  double alpha = max(alphaL,alphaR);
58  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
59  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
60  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
61  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
62  }
63  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
64  }
65 
66  return(0);
67 }
#define _Numa2DComputeSpeedofSound_(gamma, R, T0, dT, rho0, drho, EP, c)
Definition: numa2d.h:99
#define absolute(a)
Definition: math_ops.h:32
double * x
Definition: hypar.h:107
double R
Definition: numa2d.h:111
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Definition: numa2d.h:109
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
double gamma
Definition: numa2d.h:110
#define _ArrayCopy1D2_(x, y, size)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _Numa2DGetFlowVars_(u, drho, uvel, vvel, dT, rho0)
Definition: numa2d.h:33
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
#define _MODEL_NVARS_
Definition: euler1d.h:58
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
#define max(a, b)
Definition: math_ops.h:18

◆ Numa2DRusanovLinearFlux()

int Numa2DRusanovLinearFlux ( double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
double *  ,
int  ,
void *  ,
double   
)

Definition at line 69 of file Numa2DUpwind.c.

70 {
71  HyPar *solver = (HyPar*) s;
72  Numa2D *param = (Numa2D*) solver->physics;
73  int done;
74 
75  int *dim = solver->dim_local;
76  int ghosts = solver->ghosts;
77 
78  int bounds_outer[_MODEL_NDIMS_], bounds_inter[_MODEL_NDIMS_];
79  _ArrayCopy1D2_(dim,bounds_outer,_MODEL_NDIMS_); bounds_outer[dir] = 1;
80  _ArrayCopy1D2_(dim,bounds_inter,_MODEL_NDIMS_); bounds_inter[dir] += 1;
81 
82  done = 0; int index_outer[2] = {0,0}, index_inter[2];
83  while (!done) {
84  _ArrayCopy1D2_(index_outer,index_inter,_MODEL_NDIMS_);
85  for (index_inter[dir] = 0; index_inter[dir] < bounds_inter[dir]; index_inter[dir]++) {
86  int p; _ArrayIndex1D2_(_MODEL_NDIMS_,bounds_inter,index_inter,0,p);
87  double udiff[_MODEL_NVARS_],drho,vel[2],dT,rho0,P0,T0,EP,c;
88  double ycoordL, ycoordR;
89 
90  if (dir == _YDIR_) {
91  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_]-1),dim,ghosts,solver->x,ycoordL);
92  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_] ),dim,ghosts,solver->x,ycoordR);
93  } else {
94  _GetCoordinate_(_YDIR_,(index_inter[_YDIR_] ),dim,ghosts,solver->x,ycoordL);
95  ycoordR = ycoordL;
96  }
97 
98  /* Rusanov's upwinding scheme */
99 
100  udiff[0] = 0.5 * (uR[_MODEL_NVARS_*p+0] - uL[_MODEL_NVARS_*p+0]);
101  udiff[1] = 0.5 * (uR[_MODEL_NVARS_*p+1] - uL[_MODEL_NVARS_*p+1]);
102  udiff[2] = 0.5 * (uR[_MODEL_NVARS_*p+2] - uL[_MODEL_NVARS_*p+2]);
103  udiff[3] = 0.5 * (uR[_MODEL_NVARS_*p+3] - uL[_MODEL_NVARS_*p+3]);
104 
105  /* left of the interface */
106  param->StandardAtmosphere (param,ycoordL,&EP,&P0,&rho0,&T0);
107  _Numa2DGetFlowVars_ ((uL+_MODEL_NVARS_*p),drho,vel[0],vel[1],dT,rho0);
108  _Numa2DComputeLinearizedSpeedofSound_ (param->gamma,param->R,T0,rho0,EP,c);
109  double alphaL = c;
110 
111  /* right of the interface */
112  param->StandardAtmosphere (param,ycoordR,&EP,&P0,&rho0,&T0);
113  _Numa2DGetFlowVars_ ((uR+_MODEL_NVARS_*p),drho,vel[0],vel[1],dT,rho0);
114  _Numa2DComputeLinearizedSpeedofSound_ (param->gamma,param->R,T0,rho0,EP,c);
115  double alphaR = c;
116 
117  double alpha = max(alphaL,alphaR);
118  fI[_MODEL_NVARS_*p+0] = 0.5*(fL[_MODEL_NVARS_*p+0]+fR[_MODEL_NVARS_*p+0])-alpha*udiff[0];
119  fI[_MODEL_NVARS_*p+1] = 0.5*(fL[_MODEL_NVARS_*p+1]+fR[_MODEL_NVARS_*p+1])-alpha*udiff[1];
120  fI[_MODEL_NVARS_*p+2] = 0.5*(fL[_MODEL_NVARS_*p+2]+fR[_MODEL_NVARS_*p+2])-alpha*udiff[2];
121  fI[_MODEL_NVARS_*p+3] = 0.5*(fL[_MODEL_NVARS_*p+3]+fR[_MODEL_NVARS_*p+3])-alpha*udiff[3];
122 
123  /* some harmless statements to avoid compiler warnings */
124  vel[0] = vel[1];
125  dT = vel[0];
126  vel[1] = dT;
127  }
128  _ArrayIncrementIndex_(_MODEL_NDIMS_,bounds_outer,index_outer,done);
129  }
130 
131  return(0);
132 }
double * x
Definition: hypar.h:107
double R
Definition: numa2d.h:111
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
Definition: numa2d.h:109
#define _YDIR_
Definition: euler2d.h:41
#define _MODEL_NDIMS_
Definition: euler1d.h:56
double gamma
Definition: numa2d.h:110
#define _Numa2DComputeLinearizedSpeedofSound_(gamma, R, T0, rho0, EP, c)
Definition: numa2d.h:104
#define _ArrayCopy1D2_(x, y, size)
int * dim_local
Definition: hypar.h:37
#define _ArrayIncrementIndex_(N, imax, i, done)
#define _GetCoordinate_(dir, i, dim, ghosts, x, coord)
Definition: basic.h:31
void * physics
Definition: hypar.h:266
#define _Numa2DGetFlowVars_(u, drho, uvel, vvel, dT, rho0)
Definition: numa2d.h:33
int ghosts
Definition: hypar.h:52
#define _ArrayIndex1D2_(N, imax, i, ghost, index)
#define _MODEL_NVARS_
Definition: euler1d.h:58
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
#define max(a, b)
Definition: math_ops.h:18

◆ Numa2DCalculateStandardAtmosphere_1()

void Numa2DCalculateStandardAtmosphere_1 ( void *  ,
double  ,
double *  ,
double *  ,
double *  ,
double *   
)

Definition at line 5 of file Numa2DFunctions.c.

6 {
7  Numa2D *physics = (Numa2D*) p;
8 
9  double R = physics->R;
10  double gamma = physics->gamma;
11  double g = physics->g;
12 
13  /* reference quantities at zero altitude */
14  double P0, T0;
15  P0 = physics->Pref;
16  T0 = physics->Tref;
17 
18  double inv_gamma_m1 = 1.0/(gamma-1.0);
19  double Cp = gamma * inv_gamma_m1 * R;
20 
21  double theta = T0;
22  *ExnerP = 1.0 - (g/(Cp*theta))*z;
23  *P = P0 * raiseto((*ExnerP),gamma*inv_gamma_m1);
24  *rho = (P0/(R*theta)) * raiseto((*ExnerP),inv_gamma_m1);
25  *T = (*rho) * theta;
26 
27  return(0);
28 }
double Tref
Definition: numa2d.h:117
double R
Definition: numa2d.h:111
double g
Definition: numa2d.h:112
Definition: numa2d.h:109
double Pref
Definition: numa2d.h:117
double gamma
Definition: numa2d.h:110
#define raiseto(x, a)
Definition: math_ops.h:37

◆ Numa2DCalculateStandardAtmosphere_2()

void Numa2DCalculateStandardAtmosphere_2 ( void *  ,
double  ,
double *  ,
double *  ,
double *  ,
double *   
)

Definition at line 30 of file Numa2DFunctions.c.

31 {
32  Numa2D *physics = (Numa2D*) p;
33 
34  double R = physics->R;
35  double gamma = physics->gamma;
36  double g = physics->g;
37 
38  /* reference quantities at zero altitude */
39  double P0, T0;
40  P0 = physics->Pref;
41  T0 = physics->Tref;
42 
43  double BV = 0.01; /* Brunt-Vaisala frequency */
44  double inv_gamma_m1 = 1.0/(gamma-1.0);
45  double Cp = gamma * inv_gamma_m1 * R;
46 
47  double term = BV*BV*z/g;
48  double theta = T0 * exp(term);
49  *ExnerP = 1.0 + (g*g/(Cp*T0*BV*BV)) * (exp(-term) - 1.0);
50  *P = P0 * raiseto((*ExnerP),gamma*inv_gamma_m1);
51  *rho = (P0/(R*theta)) * raiseto((*ExnerP),inv_gamma_m1);
52  *T = (*rho) * theta;
53 
54  return(0);
55 }
double Tref
Definition: numa2d.h:117
double R
Definition: numa2d.h:111
double g
Definition: numa2d.h:112
Definition: numa2d.h:109
double Pref
Definition: numa2d.h:117
double gamma
Definition: numa2d.h:110
#define raiseto(x, a)
Definition: math_ops.h:37

◆ Numa2DInitialize()

int Numa2DInitialize ( void *  s,
void *  m 
)

Definition at line 25 of file Numa2DInitialize.c.

26 {
27  HyPar *solver = (HyPar*) s;
28  MPIVariables *mpi = (MPIVariables*) m;
29  Numa2D *physics = (Numa2D*) solver->physics;
30  int ferr = 0;
31 
32  static int count = 0;
33 
34  if (solver->nvars != _MODEL_NVARS_) {
35  fprintf(stderr,"Error in Numa2DInitialize(): nvars has to be %d.\n",_MODEL_NVARS_);
36  return(1);
37  }
38  if (solver->ndims != _MODEL_NDIMS_) {
39  fprintf(stderr,"Error in Numa2DInitialize(): ndims has to be %d.\n",_MODEL_NDIMS_);
40  return(1);
41  }
42 
43  /* default values */
44  physics->gamma = 1.4;
45  physics->R = 287.058; /* J kg^{-1} K^{-1} */
46  physics->g = 9.8; /* m s^{-2} */
47  physics->mu = 0.0; /* m^2 s^{-1} */
48 
49  physics->Pref = 101327.0; /* N m^{-2} */
50  physics->Tref = 288.15; /* Kelvin */
51 
52  strcpy(physics->upwind,_RUSANOV_UPWINDING_);
53 
54  /* default choice of initial atmosphere */
55  physics->init_atmos = 1;
56 
57  /* reading physical model specific inputs - rank 0 */
58  if (!mpi->rank) {
59  FILE *in;
60  if (!count) printf("Reading physical model inputs from file \"physics.inp\".\n");
61  in = fopen("physics.inp","r");
62  if (!in) printf("Warning: File \"physics.inp\" not found. Using default values.\n");
63  else {
64  char word[_MAX_STRING_SIZE_];
65  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
66  if (!strcmp(word, "begin")){
67  while (strcmp(word, "end")){
68  ferr = fscanf(in,"%s",word); if (ferr != 1) return(1);
69  if (!strcmp(word, "gamma")) {
70  ferr = fscanf(in,"%lf",&physics->gamma); if (ferr != 1) return(1);
71  } else if (!strcmp(word,"R")) {
72  ferr = fscanf(in,"%lf",&physics->R); if (ferr != 1) return(1);
73  } else if (!strcmp(word,"g")) {
74  ferr = fscanf(in,"%lf",&physics->g); if (ferr != 1) return(1);
75  } else if (!strcmp(word,"mu")) {
76  ferr = fscanf(in,"%lf",&physics->mu); if (ferr != 1) return(1);
77  } else if (!strcmp(word,"Pref")) {
78  ferr = fscanf(in,"%lf",&physics->Pref); if (ferr != 1) return(1);
79  } else if (!strcmp(word,"Tref")) {
80  ferr = fscanf(in,"%lf",&physics->Tref); if (ferr != 1) return(1);
81  } else if (!strcmp(word,"init_atmos")) {
82  ferr = fscanf(in,"%d",&physics->init_atmos); if (ferr != 1) return(1);
83  } else if (!strcmp(word,"upwinding")) {
84  ferr = fscanf(in,"%s",physics->upwind); if (ferr != 1) return(1);
85  } else if (strcmp(word,"end")) {
86  char useless[_MAX_STRING_SIZE_];
87  ferr = fscanf(in,"%s",useless); if (ferr != 1) return(ferr);
88  printf("Warning: keyword %s in file \"physics.inp\" with value %s not ",word,useless);
89  printf("recognized or extraneous. Ignoring.\n");
90  }
91  }
92  } else {
93  fprintf(stderr,"Error: Illegal format in file \"physics.inp\".\n");
94  return(1);
95  }
96  }
97  fclose(in);
98  }
99 
100 #ifndef serial
101  IERR MPIBroadcast_double (&physics->gamma ,1 ,0,&mpi->world);CHECKERR(ierr);
102  IERR MPIBroadcast_double (&physics->R ,1 ,0,&mpi->world);CHECKERR(ierr);
103  IERR MPIBroadcast_double (&physics->g ,1 ,0,&mpi->world);CHECKERR(ierr);
104  IERR MPIBroadcast_double (&physics->mu ,1 ,0,&mpi->world);CHECKERR(ierr);
105  IERR MPIBroadcast_double (&physics->Pref ,1 ,0,&mpi->world);CHECKERR(ierr);
106  IERR MPIBroadcast_double (&physics->Tref ,1 ,0,&mpi->world);CHECKERR(ierr);
107  IERR MPIBroadcast_integer (&physics->init_atmos,1 ,0,&mpi->world);CHECKERR(ierr);
109 #endif
110 
111  /* calculate the mean hydrostatic atmosphere as a function of altitude */
112  if (physics->init_atmos == 1) {
114  } else if (physics->init_atmos == 2) {
116  } else {
117  if (!mpi->rank) {
118  fprintf(stderr,"Error in Numa2DInitialize(): invalid choice of initial atmosphere (init_atmos).\n");
119  return(1);
120  }
121  }
122  CHECKERR(ierr);
123 
124  /* initializing physical model-specific functions */
125  if (!strcmp(solver->SplitHyperbolicFlux,"yes"))
126  solver->dFFunction = Numa2DStiffFlux;
127  else solver->dFFunction = NULL;
128  solver->FFunction = Numa2DFlux;
129  solver->ComputeCFL = Numa2DComputeCFL;
130  solver->SFunction = Numa2DSource;
131  if (!strcmp(physics->upwind,_RUSANOV_UPWINDING_)) {
132  solver->Upwind = Numa2DRusanovFlux;
133  if (!strcmp(solver->SplitHyperbolicFlux,"yes"))
135  else solver->UpwinddF = NULL;
136  } else {
137  if (!mpi->rank) fprintf(stderr,"Error in Numa2DInitialize(): Invalid choice of upwinding scheme.\n");
138  return(1);
139  }
140 
141  /* set the value of gamma in all the boundary objects */
142  int n;
143  DomainBoundary *boundary = (DomainBoundary*) solver->boundary;
144  for (n = 0; n < solver->nBoundaryZones; n++) boundary[n].gamma = physics->gamma;
145 
146  /* finally, hijack the main solver's dissipation function pointer
147  * to this model's own function, since it's difficult to express
148  * the dissipation terms in the general form */
150 
151  /* check that solver has the correct choice of diffusion formulation */
152  if (strcmp(solver->spatial_type_par,_NC_2STAGE_)) {
153  if (!mpi->rank) {
154  fprintf(stderr,"Error in Numa2DInitialize(): Parabolic term spatial discretization must be \"%s\"\n",_NC_2STAGE_);
155  }
156  return(1);
157  }
158 
159  count++;
160  return(0);
161 }
int Numa2DRusanovLinearFlux(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Numa2DUpwind.c:69
int init_atmos
Definition: numa2d.h:113
int nvars
Definition: hypar.h:29
void Numa2DCalculateStandardAtmosphere_2(void *, double, double *, double *, double *, double *)
double Tref
Definition: numa2d.h:117
#define IERR
Definition: basic.h:16
#define CHECKERR(ierr)
Definition: basic.h:18
int(* ParabolicFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:256
char SplitHyperbolicFlux[_MAX_STRING_SIZE_]
Definition: hypar.h:92
void * boundary
Definition: hypar.h:159
double mu
Definition: numa2d.h:114
int MPIBroadcast_integer(int *, int, int, void *)
Definition: MPIBroadcast.c:23
#define _RUSANOV_UPWINDING_
Definition: numa2d.h:130
int(* SFunction)(double *, double *, void *, void *, double)
Definition: hypar.h:317
int(* Upwind)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:295
int Numa2DParabolicFunction(double *, double *, void *, void *, double)
double R
Definition: numa2d.h:111
Structure containing the variables and function pointers defining a boundary.
int ndims
Definition: hypar.h:26
void Numa2DCalculateStandardAtmosphere_1(void *, double, double *, double *, double *, double *)
int MPIBroadcast_character(char *, int, int, void *)
Definition: MPIBroadcast.c:37
int(* FFunction)(double *, double *, int, void *, double)
Definition: hypar.h:276
int MPIBroadcast_double(double *, int, int, void *)
Definition: MPIBroadcast.c:9
double(* ComputeCFL)(void *, void *, double, double)
Definition: hypar.h:269
Structure containing all solver-specific variables and functions.
Definition: hypar.h:23
double g
Definition: numa2d.h:112
#define _MAX_STRING_SIZE_
Definition: basic.h:14
Definition: numa2d.h:109
double Pref
Definition: numa2d.h:117
#define _MODEL_NDIMS_
Definition: euler1d.h:56
double gamma
Definition: numa2d.h:110
MPI_Comm world
char spatial_type_par[_MAX_STRING_SIZE_]
Definition: hypar.h:96
int(* dFFunction)(double *, double *, int, void *, double)
Definition: hypar.h:280
void * physics
Definition: hypar.h:266
int(* UpwinddF)(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: hypar.h:300
int Numa2DStiffFlux(double *, double *, int, void *, double)
Definition: Numa2DFlux.c:43
Structure of MPI-related variables.
#define _MODEL_NVARS_
Definition: euler1d.h:58
int nBoundaryZones
Definition: hypar.h:157
void(* StandardAtmosphere)(void *, double, double *, double *, double *, double *)
Definition: numa2d.h:120
int Numa2DRusanovFlux(double *, double *, double *, double *, double *, double *, int, void *, double)
Definition: Numa2DUpwind.c:9
int Numa2DFlux(double *, double *, int, void *, double)
Definition: Numa2DFlux.c:7
int Numa2DSource(double *, double *, void *, void *, double)
Definition: Numa2DSource.c:7
char upwind[_MAX_STRING_SIZE_]
Definition: numa2d.h:123
#define _NC_2STAGE_
Definition: hypar.h:477
double Numa2DComputeCFL(void *, void *, double, double)