Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[WIP] New boundaries and floors for stability #24

Open
wants to merge 20 commits into
base: dev
Choose a base branch
from
Open

[WIP] New boundaries and floors for stability #24

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
20 commits
Select commit Hold shift + click to select a range
f3e8f53
Apply temp ceiling by lowering internal energy
Oct 26, 2020
69ed75e
Don't use ghost zones in X2 for fixup or reconstruction
Oct 26, 2020
428fb0c
Merge branch 'dev' into fix/bounds-and-floors
Oct 26, 2020
f373f4e
Merge branch 'dev' into fix/bounds-and-floors
Oct 26, 2020
1d9db2b
Turns out ghost zones exist
Nov 3, 2020
0b44dc7
Attempt at fluid-frame floors. Something in this branch broke early …
Nov 11, 2020
d9ee82e
Merge branch 'dev' of github:afd-illinois/iharm3d into fix/bounds-and…
Dec 30, 2020
eaa1529
Merge branch 'dev' of github:afd-illinois/iharm3d into fix/bounds-and…
Dec 30, 2020
a099daf
Very basic floors on post-reconstruction results.
Dec 30, 2020
08f564c
Merge branch 'fix/bounds-and-floors' of github:afd-illinois/iharm3d i…
Dec 30, 2020
31e668e
Fix trivial problem/compile stuff, add the option of complaining on b…
Dec 30, 2020
071c99b
Apply real floors after reconstruction
Dec 30, 2020
b8c5e79
The MAD problem was just an outdated torus
Dec 30, 2020
9d102dc
Consistent/example compile options for temperature ceiling applicatio…
Dec 30, 2020
1de2958
Update to Frontera default GSL on that machine
Dec 30, 2020
b9ab135
Compiling with intel 19.1.1 works on frontera, set flags for that
Dec 31, 2020
2db0843
Fix post-recon flooring. Add back PPM recon
Jan 5, 2021
ce91cdd
Make fixup callable from any zone location for appropriate post-recon…
Jan 9, 2021
3f826b3
PtoU needs a get_state call. Sigh.
Jan 9, 2021
9795ae3
When you use the state at a face, you probably wanted the state at th…
Jan 15, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 14 additions & 14 deletions core/bounds.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,11 +55,11 @@ void set_bounds(struct GridGeom *G, struct FluidState *S)
#if N1 < NG
int iactive = NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][j][iactive];
pflag[k][j][i] = pflag[k][j][iactive];
pflag[k][j][i] = 10;
#elif X1L_BOUND == OUTFLOW
int iz = 0 + NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][j][iz];
pflag[k][j][i] = pflag[k][j][iz];
pflag[k][j][i] = 10;

double rescale = G->gdet[CENT][j][iz]/G->gdet[CENT][j][i];
S->P[B1][k][j][i] *= rescale;
Expand Down Expand Up @@ -98,11 +98,11 @@ void set_bounds(struct GridGeom *G, struct FluidState *S)
#if N1 < NG
int iactive = N1 - 1 + NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][j][iactive];
pflag[k][j][i] = pflag[k][j][iactive];
pflag[k][j][i] = 10;
#elif X1R_BOUND == OUTFLOW
int iz = N1 - 1 + NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][j][iz];
pflag[k][j][i] = pflag[k][j][iz];
pflag[k][j][i] = 10;

double rescale = G->gdet[CENT][j][iz]/G->gdet[CENT][j][i];
S->P[B1][k][j][i] *= rescale;
Expand Down Expand Up @@ -147,16 +147,16 @@ void set_bounds(struct GridGeom *G, struct FluidState *S)
#if N2 < NG
int jactive = NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][jactive][i];
pflag[k][j][i] = pflag[k][jactive][i];
pflag[k][j][i] = 10;
#elif X2L_BOUND == OUTFLOW
int jz = 0 + NG ;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][jz][i];
pflag[k][j][i] = pflag[k][jz][i];
pflag[k][j][i] = 10;
#elif X2L_BOUND == POLAR
// Reflect the zone past NG by NG-j
int jrefl = NG + (NG - j) - 1;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][jrefl][i];
pflag[k][j][i] = pflag[k][jrefl][i];
pflag[k][j][i] = 10;
S->P[U2][k][j][i] *= -1.;
S->P[B2][k][j][i] *= -1.;
#endif
Expand All @@ -176,16 +176,16 @@ void set_bounds(struct GridGeom *G, struct FluidState *S)
#if N2 < NG
int jactive = N2 - 1 + NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][jactive][i];
pflag[k][j][i] = pflag[k][jactive][i];
pflag[k][j][i] = 10;
#elif X2R_BOUND == OUTFLOW
int jz = N2 - 1 + NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][jz][i];
pflag[k][j][i] = pflag[k][jz][i];
pflag[k][j][i] = 10;
#elif X2R_BOUND == POLAR
// As j grows beyond N2+NG, reflect the zone that far previous
int jrefl = (N2 + NG) + (N2 + NG - j) - 1;
PLOOP S->P[ip][k][j][i] = S->P[ip][k][jrefl][i];
pflag[k][j][i] = pflag[k][jrefl][i];
pflag[k][j][i] = 10;
S->P[U2][k][j][i] *= -1.;
S->P[B2][k][j][i] *= -1.;
#endif
Expand All @@ -209,11 +209,11 @@ void set_bounds(struct GridGeom *G, struct FluidState *S)
#if N3 < NG
int kactive = NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][kactive][j][i];
pflag[k][j][i] = pflag[kactive][j][i];
pflag[k][j][i] = 10;
#elif X3L_BOUND == OUTFLOW
int kz = 0 + NG ;
PLOOP S->P[ip][k][j][i] = S->P[ip][kz][j][i];
pflag[k][j][i] = pflag[kz][j][i];
pflag[k][j][i] = 10;
#endif
}
}
Expand All @@ -230,11 +230,11 @@ void set_bounds(struct GridGeom *G, struct FluidState *S)
#if N3 < NG
int kactive = N3-1+NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][kactive][j][i];
pflag[k][j][i] = pflag[kactive][j][i];
pflag[k][j][i] = 10;
#elif X3R_BOUND == OUTFLOW
int kz = N3 - 1 + NG;
PLOOP S->P[ip][k][j][i] = S->P[ip][kz][j][i];
pflag[k][j][i] = pflag[kz][j][i];
pflag[k][j][i] = 10;
#endif
}
}
Expand Down
6 changes: 3 additions & 3 deletions core/decs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,6 +121,7 @@
#define PPM (1)
#define WENO (2)
#define MP5 (3)
#define DONOR_CELL (4)

// Primitive and conserved variables
#define RHO (0)
Expand Down Expand Up @@ -171,15 +172,14 @@

// Failure modes
// TODO find+eliminate uses
// TODO make these the real U_to_P returns...
#define FAIL_UTOPRIM (0)
#define FAIL_VCHAR_DISCR (1)
#define FAIL_COEFF_NEG (2)
#define FAIL_COEFF_SUP (3)
#define FAIL_GAMMA (4)
#define FAIL_METRIC (5)

// U to P failure modes


// Timers
#define TIMER_RECON (1)
Expand Down Expand Up @@ -433,7 +433,7 @@ void fixup_electrons(struct FluidState *S);
#endif

// fixup.c
void fixup(struct GridGeom *G, struct FluidState *S);
void fixup(struct GridGeom *G, struct FluidState *S, int loc);
void fixup_utoprim(struct GridGeom *G, struct FluidState *S);

// fluxes.c
Expand Down
71 changes: 47 additions & 24 deletions core/fixup.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,28 +22,29 @@

static struct FluidState *Stmp;

void fixup_ceiling(struct GridGeom *G, struct FluidState *S, int i, int j, int k);
void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j, int k);
void fixup_ceiling(struct GridGeom *G, struct FluidState *S, int i, int j, int k, int loc);
void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j, int k, int loc);

// Apply floors to density, internal energy
void fixup(struct GridGeom *G, struct FluidState *S)
void fixup(struct GridGeom *G, struct FluidState *S, int loc)
{
timer_start(TIMER_FIXUP);

static int firstc = 1;
if (firstc) {Stmp = calloc(1,sizeof(struct FluidState)); firstc = 0;}

#pragma omp parallel for simd collapse(2)
ZLOOPALL fflag[k][j][i] = 0;

#pragma omp parallel for collapse(3)
ZLOOP fixup_ceiling(G, S, i, j, k);
if (firstc) {
Stmp = calloc(1,sizeof(struct FluidState));
firstc = 0;
}

// Bulk call before bsq calculation below
get_state_vec(G, S, CENT, 0, N3-1, 0, N2-1, 0, N1-1);
get_state_vec(G, S, loc, 0, N3-1, 0, N2-1, 0, N1-1);

#pragma omp parallel for collapse(3)
ZLOOP fixup_floor(G, S, i, j, k);
ZLOOP {
fflag[k][j][i] = 0;
fixup_floor(G, S, i, j, k, loc);
fixup_ceiling(G, S, i, j, k, loc);
}

// Some debug info about floors
#if DEBUG
Expand Down Expand Up @@ -86,11 +87,11 @@ void fixup(struct GridGeom *G, struct FluidState *S)
timer_stop(TIMER_FIXUP);
}

inline void fixup_ceiling(struct GridGeom *G, struct FluidState *S, int i, int j, int k)
inline void fixup_ceiling(struct GridGeom *G, struct FluidState *S, int i, int j, int k, int loc)
{
// First apply ceilings:
// 1. Limit gamma with respect to normal observer
double gamma = mhd_gamma_calc(G, S, i, j, k, CENT);
double gamma = mhd_gamma_calc(G, S, i, j, k, loc);

if (gamma > GAMMAMAX) {
fflag[k][j][i] |= HIT_FLOOR_GAMMA;
Expand All @@ -112,16 +113,27 @@ inline void fixup_ceiling(struct GridGeom *G, struct FluidState *S, int i, int j
S->P[KTOT][k][j][i] = KTOTMAX;
}
#endif

#if TEMP_ADJUST_U
// 3. Limit temperature by cooling in fluid frame
// Note this is applied after the RHO floor has taken effect
double u_ceil_temp = UORHOMAX * S->P[RHO][k][j][i];
if (S->P[UU][k][j][i] > u_ceil_temp) {
fflag[k][j][i] |= HIT_FLOOR_TEMP;
S->P[UU][k][j][i] = u_ceil_temp;
}
#endif

}

inline void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j, int k)
inline void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j, int k, int loc)
{
// Then apply floors:
// 1. Geometric hard floors, not based on fluid relationships
double rhoflr_geom, uflr_geom;
if(METRIC == MKS) {
double r, th, X[NDIM];
coord(i, j, k, CENT, X);
coord(i, j, k, loc, X);
bl_coord(X, &r, &th);

// New, steeper floor in rho
Expand All @@ -145,7 +157,7 @@ inline void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j,


// 2. Magnetic floors: impose maximum magnetization sigma = bsq/rho, inverse beta prop. to bsq/U
//get_state(G, S, i, j, k, CENT); // called above
//get_state(G, S, i, j, k, loc); // called above
double bsq = bsq_calc(S, i, j, k);
double rhoflr_b = bsq/BSQORHOMAX;
double uflr_b = bsq/BSQOUMAX;
Expand All @@ -157,6 +169,10 @@ inline void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j,
// Evaluate highest U floor
double uflr_max = MY_MAX(uflr_geom, uflr_b);

#if TEMP_ADJUST_U
// Evaluate highest RHO floor
double rhoflr_max = MY_MAX(rhoflr_geom, rhoflr_b);
#else
// 3. Temperature ceiling: impose maximum temperature
// Take floors on U into account
double rhoflr_temp = MY_MAX(S->P[UU][k][j][i] / UORHOMAX, uflr_max / UORHOMAX);
Expand All @@ -166,9 +182,14 @@ inline void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j,

// Evaluate highest RHO floor
double rhoflr_max = MY_MAX(MY_MAX(rhoflr_geom, rhoflr_b), rhoflr_temp);
#endif

if (rhoflr_max > S->P[RHO][k][j][i] || uflr_max > S->P[UU][k][j][i]) { // Apply floors

#if FLUID_FRAME_FLOORS
S->P[RHO][k][j][i] += MY_MAX(0., rhoflr_max - S->P[RHO][k][j][i]);
S->P[UU][k][j][i] += MY_MAX(0., uflr_max - S->P[UU][k][j][i]);
#else
// Initialize a dummy fluid parcel
PLOOP {
Stmp->P[ip][k][j][i] = 0;
Expand All @@ -180,12 +201,12 @@ inline void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j,
Stmp->P[UU][k][j][i] = MY_MAX(0., uflr_max - S->P[UU][k][j][i]);

// Get conserved variables for the parcel
get_state(G, Stmp, i, j, k, CENT);
prim_to_flux(G, Stmp, i, j, k, 0, CENT, Stmp->U);
get_state(G, Stmp, i, j, k, loc);
prim_to_flux(G, Stmp, i, j, k, 0, loc, Stmp->U);

// And for the current state
//get_state(G, S, i, j, k, CENT); // Called just above, or vectorized above that
prim_to_flux(G, S, i, j, k, 0, CENT, S->U);
//get_state(G, S, i, j, k, loc); // Called just above, or vectorized above that
prim_to_flux(G, S, i, j, k, 0, loc, S->U);

// Add new conserved variables to current values
PLOOP {
Expand All @@ -195,7 +216,9 @@ inline void fixup_floor(struct GridGeom *G, struct FluidState *S, int i, int j,

// Recover primitive variables
// CFG: do we get any failures here?
pflag[k][j][i] = U_to_P(G, S, i, j, k, CENT);
pflag[k][j][i] = U_to_P(G, S, i, j, k, loc);
#endif

}

#if ELECTRONS
Expand Down Expand Up @@ -283,9 +306,9 @@ void fixup_utoprim(struct GridGeom *G, struct FluidState *S)
#endif

// Make sure fixed values still abide by floors
fixup_ceiling(G, S, i, j, k);
fixup_ceiling(G, S, i, j, k, CENT);
get_state(G, S, i, j, k, CENT);
fixup_floor(G, S, i, j, k);
fixup_floor(G, S, i, j, k, CENT);
}
}

Expand Down
16 changes: 13 additions & 3 deletions core/fluxes.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,22 +119,25 @@ void lr_to_flux(struct GridGeom *G, struct FluidState *Sr,
PLOOP {
if (dir == 1) {
#pragma omp parallel for collapse(2)
ZSLOOP_REVERSE(-1, N3, -1, N2, -1, N1)
ZSLOOP_REVERSE(-1, N3, -1, N2, -1, N1) {
Sl->P[ip][k][j][i] = Sl->P[ip][k][j][i - 1];
}
} else if (dir == 2) {
#pragma omp parallel for collapse(2)
for (int k = (N3) + NG; k >= (-1) + NG; k--) {
for (int i = (N1) + NG; i >= (-1) + NG; i--) {
for (int j = (N2) + NG; j >= (-1) + NG; j--)
for (int j = (N2) + NG; j >= (-1) + NG; j--) {
Sl->P[ip][k][j][i] = Sl->P[ip][k][j - 1][i];
}
}
}
} else if (dir == 3) {
#pragma omp parallel for collapse(2)
for (int j = (N2) + NG; j >= (-1) + NG; j--) {
for (int i = (N1) + NG; i >= (-1) + NG; i--) {
for (int k = (N3) + NG; k >= (-1) + NG; k--)
for (int k = (N3) + NG; k >= (-1) + NG; k--) {
Sl->P[ip][k][j][i] = Sl->P[ip][k - 1][j][i];
}
}
}
}
Expand All @@ -144,6 +147,12 @@ void lr_to_flux(struct GridGeom *G, struct FluidState *Sr,

timer_start(TIMER_LR_STATE);

get_state_vec(G, Sl, loc, -1, N3, -1, N2, -1, N1);
get_state_vec(G, Sr, loc, -1, N3, -1, N2, -1, N1);
prim_to_flux_vec(G, Sl, 0, loc, -1, N3, -1, N2, -1, N1, Sl->U);
prim_to_flux_vec(G, Sr, 0, loc, -1, N3, -1, N2, -1, N1, Sr->U);
fixup(G, Sl, loc);
fixup(G, Sr, loc);
get_state_vec(G, Sl, loc, -1, N3, -1, N2, -1, N1);
get_state_vec(G, Sr, loc, -1, N3, -1, N2, -1, N1);

Expand All @@ -153,6 +162,7 @@ void lr_to_flux(struct GridGeom *G, struct FluidState *Sr,

timer_start(TIMER_LR_PTOF);

// TODO if we made fixup() keep lockstep, we could remove the PtoU calls
prim_to_flux_vec(G, Sl, 0, loc, -1, N3, -1, N2, -1, N1, Sl->U);
prim_to_flux_vec(G, Sl, dir, loc, -1, N3, -1, N2, -1, N1, *fluxL);

Expand Down
14 changes: 6 additions & 8 deletions core/phys.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,16 +242,16 @@ inline void mhd_vchar(struct GridGeom *G, struct FluidState *S, int i, int j, in

// Find fast magnetosonic speed
bsq = bsq_calc(S, i, j, k);
rho = fabs(S->P[RHO][k][j][i]);
u = fabs(S->P[UU][k][j][i]);
rho = S->P[RHO][k][j][i];
u = S->P[UU][k][j][i];
ef = rho + gam*u;
ee = bsq + ef;
va2 = bsq/ee;
cs2 = gam*(gam - 1.)*u/ef;

cms2 = cs2 + va2 - cs2*va2;

cms2 = (cms2 < 0) ? SMALL : cms2;
cms2 = (cms2 < SMALL) ? SMALL : cms2;
cms2 = (cms2 > 1) ? 1 : cms2;

// Require that speed of wave measured by observer q->ucon is cms2
Expand All @@ -271,15 +271,13 @@ inline void mhd_vchar(struct GridGeom *G, struct FluidState *S, int i, int j, in
B = 2.*(AuBu - (AB + AuBu)*cms2);
C = Au2 - (Asq + Au2)*cms2;

discr = B*B - 4.*A*C;
discr = (discr < 0.) ? 0. : discr;
discr = sqrt(discr);
discr = sqrt(MY_MAX(B*B - 4.*A*C, 0.));

vp = -(-B + discr)/(2.*A);
vm = -(-B - discr)/(2.*A);

cmax[k][j][i] = (vp > vm) ? vp : vm;
cmin[k][j][i] = (vp > vm) ? vm : vp;
cmax[k][j][i] = MY_MAX(vp, vm);
cmin[k][j][i] = MY_MIN(vp, vm);
}

// Source terms for equations of motion
Expand Down
Loading