In
|
if (grid_geometry == 3) then ! Cylindrical Coordinates |
|
$:GPU_PARALLEL_LOOP(collapse=4,private='[j, k, l, q, inv_ds, velocity_val, flux_face1, flux_face2]') |
|
do j = 1, sys_size |
|
do k = 0, p |
|
do q = 0, n |
|
do l = 0, m |
|
inv_ds = 1._wp/(dz(k)*y_cc(q)) |
|
velocity_val = q_prim_vf%vf(eqn_idx%cont%end + idir)%sf(l, q, k) |
|
flux_face1 = flux_n(3)%vf(j)%sf(l, q, k - 1) |
|
flux_face2 = flux_n(3)%vf(j)%sf(l, q, k) |
|
rhs_vf(j)%sf(l, q, k) = rhs_vf(j)%sf(l, q, k) + inv_ds*velocity_val*(flux_face1 - flux_face2) |
|
end do |
|
end do |
|
end do |
|
end do |
|
$:END_GPU_PARALLEL_LOOP() |
|
$:GPU_PARALLEL_LOOP(collapse=4,private='[j, k, l, q, flux_face1, flux_face2]') |
|
do j = 1, sys_size |
|
do k = 0, p |
|
do q = 0, n |
|
do l = 0, m |
|
flux_face1 = flux_gsrc_n(3)%vf(j)%sf(l, q, k - 1) |
|
flux_face2 = flux_gsrc_n(3)%vf(j)%sf(l, q, k) |
|
rhs_vf(j)%sf(l, q, k) = rhs_vf(j)%sf(l, q, k) - 5.e-1_wp/y_cc(q)*(flux_face1 + flux_face2) |
, angular flux divergence is multiplied by
u_theta/(r*dtheta) instead of
1/(r*dtheta). The flux already contains velocity. Angular fraction-source terms also omit
1/r. This changes angular transport and prevents uniform mixture fractions from remaining constant.
this case uses model_eqns=2, alt_soundspeed=F, and no phase relaxation, so its fraction equation is $D\alpha_i/Dt=0$. Consequently, the initially uniform fractions 0.4 and 0.6 must remain constant even when velocity divergence is nonzero. This follows from the volume fraction transport model described in Wilfong et al., section 2.1.1.
For a cylindrical cell, angular face area divided by cell volume is $1/(r\Delta\theta)$. This factor multiplies the numerical flux difference directly. The matching fraction source is $\alpha_i\partial_\theta u_\theta/r$, so both terms need the same metric to cancel for uniform fractions. An extra velocity multiplier also suppresses the angular pressure force when angular velocity is zero. Removing that multiplier and restoring the source metric, as in #1824, recovers the cylindrical conservation law.
Reproducer on CPU in FP64 with GCC 13.3:
{
"run_time_info": "T",
"m": 25,
"n": 25,
"p": 31,
"dt": 0.00025,
"t_step_start": 0,
"t_step_stop": 400,
"t_step_save": 20,
"num_patches": 1,
"model_eqns": 2,
"num_fluids": 2,
"alt_soundspeed": "F",
"mpp_lim": "F",
"mixture_err": "F",
"time_stepper": 3,
"recon_type": 1,
"weno_order": 5,
"weno_eps": 1e-16,
"mapped_weno": "F",
"null_weights": "F",
"mp_weno": "F",
"riemann_solver": 2,
"wave_speeds": 1,
"avg_state": 2,
"x_domain%beg": 0.0,
"x_domain%end": 1.0,
"bc_x%beg": -1,
"bc_x%end": -1,
"patch_icpp(1)%geometry": 10,
"patch_icpp(1)%x_centroid": 0.5,
"patch_icpp(1)%length_x": 1.0,
"patch_icpp(1)%vel(1)": 0.0,
"patch_icpp(1)%pres": 1.0,
"patch_icpp(1)%alpha(1)": 0.4,
"patch_icpp(1)%alpha_rho(1)": 0.4,
"patch_icpp(1)%alpha(2)": 0.6,
"patch_icpp(1)%alpha_rho(2)": 0.6,
"fluid_pp(1)%eos": 2,
"fluid_pp(1)%gamma": 2.5,
"fluid_pp(2)%eos": 2,
"fluid_pp(2)%gamma": 2.5,
"format": 1,
"precision": 2,
"parallel_io": "F",
"prim_vars_wrt": "T",
"cons_vars_wrt": "T",
"bubbles_euler": "F",
"bubbles_lagrange": "F",
"relax": "F",
"cyl_coord": "T",
"y_domain%beg": 0.0,
"y_domain%end": 1.0,
"z_domain%beg": 0.0,
"z_domain%end": 6.283185307179586,
"bc_y%beg": -14,
"bc_y%end": -2,
"bc_z%beg": -1,
"bc_z%end": -1,
"patch_icpp(1)%y_centroid": 0.0,
"patch_icpp(1)%z_centroid": 0.0,
"patch_icpp(1)%radius": 1.0,
"patch_icpp(1)%length_y": -1000000.0,
"patch_icpp(1)%length_z": -1000000.0,
"patch_icpp(1)%vel(2)": 0.0,
"patch_icpp(1)%vel(3)": "3.0*y*(1.0+0.25*sin(z))"
}./mfc.sh run case.json -t pre_process simulation --no-mpi --no-gpu -n 1 -j 8
Check the last columns of D/cons.7.00.000400.dat and D/cons.8.00.000400.dat, which contain the two volume fractions at time 0.1. They should remain 0.4 and 0.6 throughout the domain. The original code produces fluid 1 fractions from 0.3555261512299 to 0.47866057440201. With #1824, every saved fluid 1 fraction remains 0.4.
The supporting HLL rotation case checks the momentum geometry added in #1824. Initially, pressure is uniform, radial velocity is zero and tangential velocity equals radius. The initial outward acceleration should equal radius.
{
"run_time_info": "T",
"m": 25,
"n": 25,
"p": 31,
"dt": 0.0001,
"t_step_start": 0,
"t_step_stop": 1,
"t_step_save": 1,
"num_patches": 1,
"model_eqns": 2,
"num_fluids": 2,
"alt_soundspeed": "F",
"mpp_lim": "F",
"mixture_err": "F",
"time_stepper": 1,
"recon_type": 1,
"weno_order": 5,
"weno_eps": 1e-16,
"mapped_weno": "F",
"null_weights": "F",
"mp_weno": "F",
"riemann_solver": 1,
"wave_speeds": 1,
"avg_state": 2,
"x_domain%beg": 0.0,
"x_domain%end": 1.0,
"bc_x%beg": -1,
"bc_x%end": -1,
"patch_icpp(1)%geometry": 10,
"patch_icpp(1)%x_centroid": 0.5,
"patch_icpp(1)%length_x": 1.0,
"patch_icpp(1)%vel(1)": 0.0,
"patch_icpp(1)%pres": 1.0,
"patch_icpp(1)%alpha(1)": 0.4,
"patch_icpp(1)%alpha_rho(1)": 0.4,
"patch_icpp(1)%alpha(2)": 0.6,
"patch_icpp(1)%alpha_rho(2)": 0.6,
"fluid_pp(1)%eos": 2,
"fluid_pp(1)%gamma": 2.5,
"fluid_pp(2)%eos": 2,
"fluid_pp(2)%gamma": 2.5,
"format": 1,
"precision": 2,
"parallel_io": "F",
"prim_vars_wrt": "T",
"cons_vars_wrt": "T",
"bubbles_euler": "F",
"bubbles_lagrange": "F",
"relax": "F",
"cyl_coord": "T",
"y_domain%beg": 0.0,
"y_domain%end": 1.0,
"z_domain%beg": 0.0,
"z_domain%end": 6.283185307179586,
"bc_y%beg": -14,
"bc_y%end": -2,
"bc_z%beg": -1,
"bc_z%end": -1,
"patch_icpp(1)%y_centroid": 0.0,
"patch_icpp(1)%z_centroid": 0.0,
"patch_icpp(1)%radius": 1.0,
"patch_icpp(1)%length_y": -1000000.0,
"patch_icpp(1)%length_z": -1000000.0,
"patch_icpp(1)%vel(2)": 0.0,
"patch_icpp(1)%vel(3)": "y"
}./mfc.sh run case.json -t pre_process simulation --no-mpi --no-gpu -n 1 -j 8
In the output at step 1, compute radial velocity as cons.4 / (cons.1+cons.2) and divide by the time step 0.0001. Compare with the radial cell center, obtained by subtracting 1/52 from the second coordinate column. The original code gives zero outward acceleration. With #1824, the maximum absolute difference from the expected radius is below 0.000000000047.
Golden case IDs: 301B9153, 128954AD, 07C33719, 939D6718, 09623DE3.
In
MFC/src/simulation/m_rhs.fpp
Lines 1312 to 1335 in 511cda5
u_theta/(r*dtheta)instead of1/(r*dtheta). The flux already contains velocity. Angular fraction-source terms also omit1/r. This changes angular transport and prevents uniform mixture fractions from remaining constant.this case uses$D\alpha_i/Dt=0$ . Consequently, the initially uniform fractions 0.4 and 0.6 must remain constant even when velocity divergence is nonzero. This follows from the volume fraction transport model described in Wilfong et al., section 2.1.1.
model_eqns=2,alt_soundspeed=F, and no phase relaxation, so its fraction equation isFor a cylindrical cell, angular face area divided by cell volume is$1/(r\Delta\theta)$ . This factor multiplies the numerical flux difference directly. The matching fraction source is $\alpha_i\partial_\theta u_\theta/r$ , so both terms need the same metric to cancel for uniform fractions. An extra velocity multiplier also suppresses the angular pressure force when angular velocity is zero. Removing that multiplier and restoring the source metric, as in #1824, recovers the cylindrical conservation law.
Reproducer on CPU in FP64 with GCC 13.3:
{ "run_time_info": "T", "m": 25, "n": 25, "p": 31, "dt": 0.00025, "t_step_start": 0, "t_step_stop": 400, "t_step_save": 20, "num_patches": 1, "model_eqns": 2, "num_fluids": 2, "alt_soundspeed": "F", "mpp_lim": "F", "mixture_err": "F", "time_stepper": 3, "recon_type": 1, "weno_order": 5, "weno_eps": 1e-16, "mapped_weno": "F", "null_weights": "F", "mp_weno": "F", "riemann_solver": 2, "wave_speeds": 1, "avg_state": 2, "x_domain%beg": 0.0, "x_domain%end": 1.0, "bc_x%beg": -1, "bc_x%end": -1, "patch_icpp(1)%geometry": 10, "patch_icpp(1)%x_centroid": 0.5, "patch_icpp(1)%length_x": 1.0, "patch_icpp(1)%vel(1)": 0.0, "patch_icpp(1)%pres": 1.0, "patch_icpp(1)%alpha(1)": 0.4, "patch_icpp(1)%alpha_rho(1)": 0.4, "patch_icpp(1)%alpha(2)": 0.6, "patch_icpp(1)%alpha_rho(2)": 0.6, "fluid_pp(1)%eos": 2, "fluid_pp(1)%gamma": 2.5, "fluid_pp(2)%eos": 2, "fluid_pp(2)%gamma": 2.5, "format": 1, "precision": 2, "parallel_io": "F", "prim_vars_wrt": "T", "cons_vars_wrt": "T", "bubbles_euler": "F", "bubbles_lagrange": "F", "relax": "F", "cyl_coord": "T", "y_domain%beg": 0.0, "y_domain%end": 1.0, "z_domain%beg": 0.0, "z_domain%end": 6.283185307179586, "bc_y%beg": -14, "bc_y%end": -2, "bc_z%beg": -1, "bc_z%end": -1, "patch_icpp(1)%y_centroid": 0.0, "patch_icpp(1)%z_centroid": 0.0, "patch_icpp(1)%radius": 1.0, "patch_icpp(1)%length_y": -1000000.0, "patch_icpp(1)%length_z": -1000000.0, "patch_icpp(1)%vel(2)": 0.0, "patch_icpp(1)%vel(3)": "3.0*y*(1.0+0.25*sin(z))" }Check the last columns of
D/cons.7.00.000400.datandD/cons.8.00.000400.dat, which contain the two volume fractions at time 0.1. They should remain 0.4 and 0.6 throughout the domain. The original code produces fluid 1 fractions from 0.3555261512299 to 0.47866057440201. With #1824, every saved fluid 1 fraction remains 0.4.The supporting HLL rotation case checks the momentum geometry added in #1824. Initially, pressure is uniform, radial velocity is zero and tangential velocity equals radius. The initial outward acceleration should equal radius.
{ "run_time_info": "T", "m": 25, "n": 25, "p": 31, "dt": 0.0001, "t_step_start": 0, "t_step_stop": 1, "t_step_save": 1, "num_patches": 1, "model_eqns": 2, "num_fluids": 2, "alt_soundspeed": "F", "mpp_lim": "F", "mixture_err": "F", "time_stepper": 1, "recon_type": 1, "weno_order": 5, "weno_eps": 1e-16, "mapped_weno": "F", "null_weights": "F", "mp_weno": "F", "riemann_solver": 1, "wave_speeds": 1, "avg_state": 2, "x_domain%beg": 0.0, "x_domain%end": 1.0, "bc_x%beg": -1, "bc_x%end": -1, "patch_icpp(1)%geometry": 10, "patch_icpp(1)%x_centroid": 0.5, "patch_icpp(1)%length_x": 1.0, "patch_icpp(1)%vel(1)": 0.0, "patch_icpp(1)%pres": 1.0, "patch_icpp(1)%alpha(1)": 0.4, "patch_icpp(1)%alpha_rho(1)": 0.4, "patch_icpp(1)%alpha(2)": 0.6, "patch_icpp(1)%alpha_rho(2)": 0.6, "fluid_pp(1)%eos": 2, "fluid_pp(1)%gamma": 2.5, "fluid_pp(2)%eos": 2, "fluid_pp(2)%gamma": 2.5, "format": 1, "precision": 2, "parallel_io": "F", "prim_vars_wrt": "T", "cons_vars_wrt": "T", "bubbles_euler": "F", "bubbles_lagrange": "F", "relax": "F", "cyl_coord": "T", "y_domain%beg": 0.0, "y_domain%end": 1.0, "z_domain%beg": 0.0, "z_domain%end": 6.283185307179586, "bc_y%beg": -14, "bc_y%end": -2, "bc_z%beg": -1, "bc_z%end": -1, "patch_icpp(1)%y_centroid": 0.0, "patch_icpp(1)%z_centroid": 0.0, "patch_icpp(1)%radius": 1.0, "patch_icpp(1)%length_y": -1000000.0, "patch_icpp(1)%length_z": -1000000.0, "patch_icpp(1)%vel(2)": 0.0, "patch_icpp(1)%vel(3)": "y" }In the output at step 1, compute radial velocity as
cons.4 / (cons.1+cons.2)and divide by the time step 0.0001. Compare with the radial cell center, obtained by subtracting 1/52 from the second coordinate column. The original code gives zero outward acceleration. With #1824, the maximum absolute difference from the expected radius is below 0.000000000047.Golden case IDs:
301B9153,128954AD,07C33719,939D6718,09623DE3.