-
Notifications
You must be signed in to change notification settings - Fork 1
/
vhtimpl.h
179 lines (169 loc) · 6.4 KB
/
vhtimpl.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
#if !defined _vhtcase_h
#define _vhtcase_h
#include <dohpfs.h>
#include <dohpvec.h>
#include <dohpsys.h>
#include <dohpunits.h>
#include <dohp.h>
struct VHTRheology {
dReal B0; /* Viscosity at reference strain rate and temperature */
dReal Bomega; /* Softening due to water content */
dReal R; /* Ideal gas constant, for Arrhenius relation */
dReal Q; /* "Activation energy" for creep */
dReal V; /* "Activation volume" for creep */
dReal du0; /* Reference strain rate */
dReal gamma0; /* Second invariant of reference strain rate */
dReal eps; /* Fraction of reference strain rate at which to regularize */
dReal pe; /* Rheological exponent */
dReal k_T; /* Thermal conductivity in the cold part */
dReal kappa_w; /* Hydraulic conductivity in the warm part */
dReal c_i; /* Specific heat capacity of cold part */
dReal Latent; /* Latent heat of fusion */
dReal rhoi; /* Density of cold part */
dReal rhow; /* Density of melt */
dReal beta_CC; /* Clausius-Capeyron gradient */
dReal T0; /* Reference temperature for definition of enthalpy and viscosity */
dReal p0; /* Reference pressure, model pressure is the deviation from reference pressure */
dReal T3; /* Triple point temperature */
dReal splice_delta; /* Characteristic width of splice */
dReal gravity[3]; /* Gravity vector */
dReal Kstab; /* Stabilization for energy diffusion */
dReal supg; /* Multiplier for SU/PG stabilization */
dReal supg_crosswind; /* Fraction of streamline diffusion to put in the cross-wind direction */
dReal shockhmm; /* Multiplier for Hughes-Mallet-Mizukami shock-capturing term */
dReal expstab; /* Multiplier for the exponential low-temperature stabilization */
dReal expstab_rate; /* Inverse characteristic width of exponential low-temperature stabilization, units of 1./(T3-T0) */
dReal mask_kinetic; /* Parameter to turn on the use of kinetic energy when computing velocity */
dReal mask_momtrans; /* Multiplier for the transport term in momentum balance */
dReal mask_rho; /* Multiplier for the true rho */
dReal mask_Ep; /* Multiplier for p in (E+p) term in energy equation */
dReal eta_min,eta_max;
dReal small; /* Small number used to avoid dividing by zero */
struct {
dReal momentum,mass,energy;
} rscale; /* Residual scaling for momentum,mass,energy */
};
struct VHTUnitTable {
dUnit Length;
dUnit Mass;
dUnit Time;
dUnit Temperature;
dUnit Density;
dUnit Energy;
dUnit Pressure;
dUnit StrainRate;
dUnit Velocity;
dUnit Acceleration;
dUnit Viscosity;
dUnit Volume;
dUnit EnergyDensity;
dUnit MomentumDensity;
dUnit EnergyPerTemperature;
dUnit ThermalConductivity;
dUnit HydroConductivity;
dUnit Diffusivity;
dUnit SpecificHeat;
dUnit EnergyPerMass;
dUnit CCGradient;
};
typedef struct _n_VHTCase *VHTCase;
struct _n_VHTCase {
MPI_Comm comm;
dErr (*solution)(VHTCase,const dReal x[3],dScalar u[],dScalar du[],dScalar *p,dScalar dp[],dScalar *e,dScalar *de);
dErr (*forcing)(VHTCase,const dReal x[3],dScalar fu[],dScalar *fp,dScalar *fe);
dErr (*setfromoptions)(VHTCase);
dErr (*destroy)(VHTCase);
struct VHTRheology rheo;
dReal bbox[3][2];
dBool reality; // The "solution" is just a guess or boundary conditions
dUnits units;
struct VHTUnitTable utable;
char name[dNAME_LEN];
void *data;
};
typedef dErr (*VHTCaseCreateFunction)(VHTCase);
struct VHTPack {
dScalar rhou[3],p,E;
dScalar drhou[9],dp[3],dE[3];
};
typedef struct {
dReal x;
dReal dp,dE;
} VHTScalarD;
struct VHTStash {
dReal Dui[6]; /* Symmetrized velocity gradient */
dReal u[3]; /* Total Velocity */
dReal K[2]; /* Total pressure and energy diffusivities, -K[0]*dp - K[1]*dE */
dReal K1[2][2]; /* Derivatives */
VHTScalarD T;
VHTScalarD omega;
dScalar omega2[3]; /* Second derivatives: o2pp, o2pE=o2Ep, o2EE */
VHTScalarD eta;
dReal eta1gamma;
VHTScalarD rho;
dReal p,E;
dReal dE[3];
dReal dp[3];
};
struct VHTLogEpoch {
dReal eta[2];
dReal cPeclet[2];
dReal cReynolds[2];
dReal p[2];
dReal E[2];
dReal T[2];
dReal K1[2];
dReal omega[2];
dReal Prandtl[2];
dReal upwind[2];
dReal speed[2];
};
struct VHTLog {
dInt epoch,alloc;
dBool monitor;
struct VHTLogEpoch *epochs;
struct VHTLogEpoch global;
};
typedef enum {EVAL_FUNCTION,EVAL_JACOBIAN, EVAL_UB} VHTEvaluation;
typedef enum {VHT_MULT_UU,VHT_MULT_UP,VHT_MULT_UE,VHT_MULT_PU,VHT_MULT_EU,VHT_MULT_EP,VHT_MULT_EE} VHTMultMode;
typedef struct _n_VHT *VHT;
struct _n_VHT {
MPI_Comm comm;
VHTCase scase;
dFS fsu,fsp,fse;
Vec xu,xp,xe,yu,yp,ye;
Vec gvelocity,gpressure,genergy;
Vec gpacked;
struct {
IS ublock,pblock;
IS lublock,lpblock;
VecScatter extractVelocity,extractPressure;
Vec x,y;
} stokes;
struct {
IS ublock,pblock,eblock; /* Global index sets for each block */
IS lublock,lpblock,leblock; /* Local index sets for each block */
IS sblock,lsblock;
VecScatter extractVelocity,extractPressure,extractEnergy;
} all;
dInt velocityBDeg,pressureCodim,energyBDeg;
dBool cardinalMass;
char mattype_Buu[256],mattype_Bpp[256],mattype_Bee[256];
dQuadratureMethod function_qmethod,jacobian_qmethod;
dRulesetIterator regioniter[EVAL_UB];
dInt u_dirichlet[16]; /* Set numbers for Dirichlet conditions, 0 means unused */
dInt e_dirichlet[16]; /* Set numbers for Dirichlet conditions, 0 means unused */
dBool alldirichlet;
dBool split_recursive;
dInt domain_error;
dBool function_monitor;
struct {
dReal p[2],E[2];
} clip;
struct VHTLog log;
};
extern PetscFunctionList VHTCaseList;
dErr VHTCaseRegister(const char *name,VHTCaseCreateFunction screate);
dErr VHTCaseRegisterAll_Exact(void); /* Defined in generated vhtexact.c */
dErr VHTCaseRegisterAll_Jako(void);
#endif