Skip to content

Commit

Permalink
Merge branch 'dev' into dev-mhdSlicesOutput
Browse files Browse the repository at this point in the history
  • Loading branch information
@evaneschneider
evaneschneider authored Sep 14, 2023
2 parents 9ab72fe + e8021cb commit 6609d9e
Show file tree
Hide file tree
Showing 2 changed files with 71 additions and 65 deletions.
2 changes: 1 addition & 1 deletion examples/3D/KH_res_ind_3D.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ ny=128
# number of grid cells in the z dimension
nz=128
# final output time
tout=5.0
tout=3.0
# time interval for output
outstep=0.01
# value of gamma
Expand Down
134 changes: 70 additions & 64 deletions src/grid/initial_conditions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -772,8 +772,8 @@ void Grid3D::KH_res_ind()

d1 = 100.0; // inner density
d2 = 1.0; // outer density
v1 = 10.5; // inner velocity
v2 = 9.5; // outer velocity
v1 = 0.5; // inner velocity
v2 = -0.5; // outer velocity
P = 2.5; // pressure
dy = 0.05; // width of ramp function (see Robertson 2009)
A = 0.1; // amplitude of the perturbation
Expand All @@ -788,73 +788,79 @@ void Grid3D::KH_res_ind()
id = i + j * H.nx + k * H.nx * H.ny;
// get the centered x and y positions
Get_Position(i, j, k, &x_pos, &y_pos, &z_pos);

// inner fluid
if (fabs(y_pos - 0.5) < 0.25) {
if (y_pos > 0.5) {
C.density[id] =
d1 - (d1 - d2) * exp(-0.5 * pow(y_pos - 0.75 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] =
v1 * C.density[id] - C.density[id] * (v1 - v2) *
exp(-0.5 * pow(y_pos - 0.75 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(y_pos - 0.75 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
} else {
C.density[id] =
d1 - (d1 - d2) * exp(-0.5 * pow(y_pos - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] =
v1 * C.density[id] - C.density[id] * (v1 - v2) *
exp(-0.5 * pow(y_pos - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(y_pos - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
// 2D initial conditions:
if (H.nz == 1) {
// inner fluid
if (fabs(y_pos - 0.5) < 0.25) {
if (y_pos > 0.5) {
C.density[id] =
d1 - (d1 - d2) * exp(-0.5 * pow(y_pos - 0.75 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] = v1 * C.density[id] -
C.density[id] * (v1 - v2) *
exp(-0.5 * pow(y_pos - 0.75 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(y_pos - 0.75 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
} else {
C.density[id] =
d1 - (d1 - d2) * exp(-0.5 * pow(y_pos - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] = v1 * C.density[id] -
C.density[id] * (v1 - v2) *
exp(-0.5 * pow(y_pos - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(y_pos - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
}
}
}
// outer fluid
else {
if (y_pos > 0.5) {
C.density[id] =
d2 + (d1 - d2) * exp(-0.5 * pow(y_pos - 0.75 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
// outer fluid
else {
if (y_pos > 0.5) {
C.density[id] =
d2 + (d1 - d2) * exp(-0.5 * pow(y_pos - 0.75 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] = v2 * C.density[id] +
C.density[id] * (v1 - v2) *
exp(-0.5 * pow(y_pos - 0.75 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(y_pos - 0.75 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
} else {
C.density[id] =
d2 + (d1 - d2) * exp(-0.5 * pow(y_pos - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] = v2 * C.density[id] +
C.density[id] * (v1 - v2) *
exp(-0.5 * pow(y_pos - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(y_pos - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
}
}
// C.momentum_y[id] = C.density[id] * A*sin(4*PI*x_pos);
C.momentum_z[id] = 0.0;

// 3D initial conditions:
} else {
// cylindrical version (3D only)
r = sqrt((z_pos - zc) * (z_pos - zc) + (y_pos - yc) * (y_pos - yc)); // center the cylinder at yc, zc
phi = atan2((z_pos - zc), (y_pos - yc));

if (r < 0.25) // inside the cylinder
{
C.density[id] = d1 - (d1 - d2) * exp(-0.5 * pow(r - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] =
v2 * C.density[id] + C.density[id] * (v1 - v2) *
exp(-0.5 * pow(y_pos - 0.75 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(y_pos - 0.75 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
} else {
C.density[id] =
d2 + (d1 - d2) * exp(-0.5 * pow(y_pos - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
v1 * C.density[id] -
C.density[id] * exp(-0.5 * pow(r - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = cos(phi) * C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_z[id] = sin(phi) * C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
} else // outside the cylinder
{
C.density[id] = d2 + (d1 - d2) * exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] =
v2 * C.density[id] + C.density[id] * (v1 - v2) *
exp(-0.5 * pow(y_pos - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(y_pos - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
v2 * C.density[id] +
C.density[id] * exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = cos(phi) * C.density[id] * A * sin(4 * M_PI * x_pos) *
(1.0 - exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy)));
C.momentum_z[id] = sin(phi) * C.density[id] * A * sin(4 * M_PI * x_pos) *
(1.0 - exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy)));
}
}
// C.momentum_y[id] = C.density[id] * A*sin(4*PI*x_pos);
C.momentum_z[id] = 0.0;

// cylindrical version (3D only)
r = sqrt((z_pos - zc) * (z_pos - zc) + (y_pos - yc) * (y_pos - yc)); // center the cylinder at yc, zc
phi = atan2((z_pos - zc), (y_pos - yc));

if (r < 0.25) // inside the cylinder
{
C.density[id] = d1 - (d1 - d2) * exp(-0.5 * pow(r - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] = v1 * C.density[id] -
C.density[id] * exp(-0.5 * pow(r - 0.25 - sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = cos(phi) * C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_z[id] = sin(phi) * C.density[id] * A * sin(4 * M_PI * x_pos) *
exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
} else // outside the cylinder
{
C.density[id] = d2 + (d1 - d2) * exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_x[id] = v2 * C.density[id] +
C.density[id] * exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy));
C.momentum_y[id] = cos(phi) * C.density[id] * A * sin(4 * M_PI * x_pos) *
(1.0 - exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy)));
C.momentum_z[id] = sin(phi) * C.density[id] * A * sin(4 * M_PI * x_pos) *
(1.0 - exp(-0.5 * pow(r - 0.25 + sqrt(-2.0 * dy * dy * log(0.5)), 2) / (dy * dy)));
}

// No matter what we do with the density and momentum, set the Energy
// and GasEnergy appropriately
Expand Down

0 comments on commit 6609d9e

Please sign in to comment.