Bug: cylindrical angular flux has an extra velocity factor and missing metric scaling

[aliencaocao]

In

MFC/src/simulation/m_rhs.fpp

Lines 1312 to 1335 in 511cda5

	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.

Minimal repro:
Based on 301B9153 but use the case.json below:

{
  "run_time_info": "T",
  "m": 29,
  "n": 29,
  "p": 29,
  "dt": 0.0005,
  "t_step_start": 0,
  "t_step_stop": 1,
  "t_step_save": 1,
  "num_patches": 3,
  "model_eqns": 2,
  "alt_soundspeed": "F",
  "num_fluids": 2,
  "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,
  "format": 1,
  "precision": 2,
  "patch_icpp(1)%pres": 1.0,
  "patch_icpp(1)%alpha_rho(1)": 0.81,
  "patch_icpp(1)%alpha(1)": 0.4,
  "patch_icpp(2)%pres": 0.5,
  "patch_icpp(2)%alpha_rho(1)": 0.25,
  "patch_icpp(2)%alpha(1)": 0.4,
  "patch_icpp(3)%pres": 0.1,
  "patch_icpp(3)%alpha_rho(1)": 0.08,
  "patch_icpp(3)%alpha(1)": 0.4,
  "fluid_pp(1)%gamma": 2.5000000000000004,
  "fluid_pp(1)%eos": 2,
  "fluid_pp(1)%cv": 0.0,
  "fluid_pp(1)%qv": 0.0,
  "fluid_pp(1)%qvp": 0.0,
  "bubbles_euler": "F",
  "bubble_model": 3,
  "polytropic": "T",
  "polydisperse": "F",
  "thermal": 3,
  "patch_icpp(1)%r0": 1,
  "patch_icpp(1)%v0": 0,
  "patch_icpp(2)%r0": 1,
  "patch_icpp(2)%v0": 0,
  "patch_icpp(3)%r0": 1,
  "patch_icpp(3)%v0": 0,
  "qbmm": "F",
  "dist_type": 2,
  "poly_sigma": 0.3,
  "sigR": 0.1,
  "sigV": 0.1,
  "rhoRV": 0.0,
  "acoustic_source": "F",
  "num_source": 1,
  "acoustic(1)%loc(1)": 0.5,
  "acoustic(1)%mag": 0.2,
  "acoustic(1)%length": 0.25,
  "acoustic(1)%dir": 1.0,
  "acoustic(1)%npulse": 1,
  "acoustic(1)%pulse": 1,
  "rdma_mpi": "F",
  "bubbles_lagrange": "F",
  "lag_params%nBubs_glb": 1,
  "lag_params%solver_approach": 0,
  "lag_params%cluster_type": 2,
  "lag_params%pressure_corrector": "F",
  "lag_params%smooth_type": 1,
  "lag_params%epsilonb": 1.0,
  "lag_params%heatTransfer_model": "F",
  "lag_params%massTransfer_model": "F",
  "lag_params%valmaxvoid": 0.9,
  "x_domain%beg": 0.0,
  "x_domain%end": 5.0,
  "y_domain%beg": 0.0,
  "y_domain%end": 1.0,
  "z_domain%beg": 0.0,
  "z_domain%end": 6.283185307179586,
  "bc_x%beg": -3,
  "bc_x%end": -3,
  "bc_y%beg": -14,
  "bc_y%end": -3,
  "bc_z%beg": -1,
  "bc_z%end": -1,
  "patch_icpp(1)%geometry": 10,
  "patch_icpp(1)%z_centroid": 0.0,
  "patch_icpp(1)%length_z": -1000000.0,
  "patch_icpp(2)%z_centroid": 0.0,
  "patch_icpp(2)%length_z": -1000000.0,
  "patch_icpp(3)%z_centroid": 0.0,
  "patch_icpp(3)%length_z": -1000000.0,
  "patch_icpp(1)%y_centroid": 0.0,
  "patch_icpp(1)%length_y": -1000000.0,
  "patch_icpp(1)%x_centroid": 0.5,
  "patch_icpp(1)%length_x": 1.0,
  "patch_icpp(1)%vel(1)": 0.0,
  "patch_icpp(1)%vel(2)": "0.1*y",
  "patch_icpp(1)%vel(3)": "0.2*y*(1.0+0.1*sin(z))",
  "patch_icpp(2)%geometry": 10,
  "patch_icpp(2)%y_centroid": 0.0,
  "patch_icpp(2)%length_y": -1000000.0,
  "patch_icpp(2)%x_centroid": 2.5,
  "patch_icpp(2)%length_x": 3.0,
  "patch_icpp(2)%vel(1)": 0.0,
  "patch_icpp(2)%vel(2)": "0.1*y",
  "patch_icpp(2)%vel(3)": "0.2*y*(1.0+0.1*sin(z))",
  "patch_icpp(3)%geometry": 10,
  "patch_icpp(3)%y_centroid": 0.0,
  "patch_icpp(3)%length_y": -1000000.0,
  "patch_icpp(3)%x_centroid": 4.5,
  "patch_icpp(3)%length_x": 1.0,
  "patch_icpp(3)%vel(1)": 0.0,
  "patch_icpp(3)%vel(2)": "0.1*y",
  "patch_icpp(3)%vel(3)": "0.2*y*(1.0+0.1*sin(z))",
  "cyl_coord": "T",
  "patch_icpp(1)%radius": 1.0,
  "patch_icpp(2)%radius": 1.0,
  "patch_icpp(3)%radius": 1.0,
  "fluid_pp(2)%gamma": 2.5,
  "fluid_pp(2)%eos": 2,
  "patch_icpp(1)%alpha_rho(2)": 0.19,
  "patch_icpp(1)%alpha(2)": 0.6,
  "patch_icpp(2)%alpha_rho(2)": 0.25,
  "patch_icpp(2)%alpha(2)": 0.6,
  "patch_icpp(3)%alpha_rho(2)": 0.0225,
  "patch_icpp(3)%alpha(2)": 0.6,
  "parallel_io": "F",
  "prim_vars_wrt": "T",
  "cons_vars_wrt": "T"
}

./mfc.sh build -t pre_process simulation --no-mpi --no-gpu -j "$(nproc)"
./mfc.sh run cylindrical/case.json -t pre_process simulation --no-mpi --no-gpu -n 1 -j "$(nproc)"

Check the last columns of cylindrical/D/cons.7.00.000001.dat and cons.8.00.000001.dat. They should remain 0.4 and 0.6. Now, it changes them by up to 4.7062e-6 after one step.

I have a fix available but I first would like to get @wilfonba 's confirmation that this is indeed a bug, because there are no golden cases that exercises this. The original cylindrical test has no initial angular velocity or angular variation.

Fix would be to remove the extra velocity and apply 1/(r*dtheta) consistently.

[wilfonba]

@aliencaocao I have some other cylindrical tests that have been giving me issues for a while now, so I'll make this change and see if it fixes those cases.

[aliencaocao]

I have a fix ready, if you would like to save some time i can PR directly.

[wilfonba]

Feel free to open the PR as a draft, and I'll run some tests over the next week or so to see whether it resolves the problem I was dealing with before.