From 4ea78a30ad11a5fafeb0ec0d5f7161256c2e1e40 Mon Sep 17 00:00:00 2001 From: ezlimen Date: Fri, 8 Sep 2023 11:34:47 -0700 Subject: [PATCH 1/4] fix builds --- builds/make.host.lux | 2 +- builds/setup.lux.sh | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/builds/make.host.lux b/builds/make.host.lux index edf4e42c0..6ce455fcb 100644 --- a/builds/make.host.lux +++ b/builds/make.host.lux @@ -9,7 +9,7 @@ GPUFLAGS = -std=c++17 OMP_NUM_THREADS = 10 #-- Library -CUDA_ROOT = /cm/shared/apps/cuda10.2/toolkit/current +CUDA_ROOT = /cm/shared/apps/cuda11.2/toolkit/current HDF5_ROOT = /cm/shared/apps/hdf5/1.10.6 FFTW_ROOT = /home/brvillas/code/fftw-3.3.8 PFFT_ROOT = /data/groups/comp-astro/bruno/code_mpi_local/pfft diff --git a/builds/setup.lux.sh b/builds/setup.lux.sh index 6d6d408f3..3ef07c50c 100755 --- a/builds/setup.lux.sh +++ b/builds/setup.lux.sh @@ -1,6 +1,7 @@ #!/bin/bash -module load hdf5/1.10.6 cuda10.2/10.2 openmpi/4.0.1 +###module load hdf5/1.10.6 cuda10.2/10.2 openmpi/4.0.1 +module load hdf5/1.10.6 cuda11.2 openmpi/4.0.1 devtoolset-9 export MACHINE=lux export CHOLLA_ENVSET=1 From 1a3ae6667bdcf8650228523cf1c3528c470855c3 Mon Sep 17 00:00:00 2001 From: ezlimen Date: Sat, 9 Sep 2023 16:15:03 -0700 Subject: [PATCH 2/4] change output time of 3D KH_res_ind and make 2D kh_res_ind ICs set --- examples/3D/KH_res_ind_3D.txt | 2 +- src/grid/initial_conditions.cpp | 126 ++++++++++++++++---------------- 2 files changed, 66 insertions(+), 62 deletions(-) diff --git a/examples/3D/KH_res_ind_3D.txt b/examples/3D/KH_res_ind_3D.txt index ab846867a..2ebe6cda0 100644 --- a/examples/3D/KH_res_ind_3D.txt +++ b/examples/3D/KH_res_ind_3D.txt @@ -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 diff --git a/src/grid/initial_conditions.cpp b/src/grid/initial_conditions.cpp index 38967e4b7..eba91f463 100644 --- a/src/grid/initial_conditions.cpp +++ b/src/grid/initial_conditions.cpp @@ -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 @@ -788,73 +788,77 @@ 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] = + // 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] = + 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] = + 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] = + 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] = + 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)); - C.momentum_x[id] = + 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] = + 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] = + 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; - - // 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))); - } + 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] = 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 From d92ee47bdd98f20c3f3e71e765144d0248d19510 Mon Sep 17 00:00:00 2001 From: ezlimen Date: Mon, 11 Sep 2023 09:29:26 -0700 Subject: [PATCH 3/4] Revert "fix builds" This reverts commit 4ea78a30ad11a5fafeb0ec0d5f7161256c2e1e40. --- builds/make.host.lux | 2 +- builds/setup.lux.sh | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/builds/make.host.lux b/builds/make.host.lux index 6ce455fcb..edf4e42c0 100644 --- a/builds/make.host.lux +++ b/builds/make.host.lux @@ -9,7 +9,7 @@ GPUFLAGS = -std=c++17 OMP_NUM_THREADS = 10 #-- Library -CUDA_ROOT = /cm/shared/apps/cuda11.2/toolkit/current +CUDA_ROOT = /cm/shared/apps/cuda10.2/toolkit/current HDF5_ROOT = /cm/shared/apps/hdf5/1.10.6 FFTW_ROOT = /home/brvillas/code/fftw-3.3.8 PFFT_ROOT = /data/groups/comp-astro/bruno/code_mpi_local/pfft diff --git a/builds/setup.lux.sh b/builds/setup.lux.sh index 3ef07c50c..6d6d408f3 100755 --- a/builds/setup.lux.sh +++ b/builds/setup.lux.sh @@ -1,7 +1,6 @@ #!/bin/bash -###module load hdf5/1.10.6 cuda10.2/10.2 openmpi/4.0.1 -module load hdf5/1.10.6 cuda11.2 openmpi/4.0.1 devtoolset-9 +module load hdf5/1.10.6 cuda10.2/10.2 openmpi/4.0.1 export MACHINE=lux export CHOLLA_ENVSET=1 From 4063d268319ce6607a7a64f35e634d5310d310ba Mon Sep 17 00:00:00 2001 From: evazlimen <109487593+evazlimen@users.noreply.github.com> Date: Mon, 11 Sep 2023 15:36:47 -0700 Subject: [PATCH 4/4] clang-format --- src/grid/initial_conditions.cpp | 148 ++++++++++++++++---------------- 1 file changed, 75 insertions(+), 73 deletions(-) diff --git a/src/grid/initial_conditions.cpp b/src/grid/initial_conditions.cpp index eba91f463..fded9236b 100644 --- a/src/grid/initial_conditions.cpp +++ b/src/grid/initial_conditions.cpp @@ -772,8 +772,8 @@ void Grid3D::KH_res_ind() d1 = 100.0; // inner density d2 = 1.0; // outer density - v1 = 0.5; // inner velocity - v2 = -0.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 @@ -788,77 +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); - // 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)); - 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] = 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))); - } - } + // 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)); + 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] = + 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