diff --git a/case_new/2D_orszag_tang_weno3m/case.py b/case_new/2D_orszag_tang_weno3m/case.py new file mode 100644 index 0000000000..195bb5e082 --- /dev/null +++ b/case_new/2D_orszag_tang_weno3m/case.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +import json +import math + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0, + "x_domain%end": 1.0, + "y_domain%beg": 0, + "y_domain%end": 1.0, + "m": 512, + "n": 512, + "p": 0, + "dt": 0.0004, + "t_step_start": 0, + "t_step_stop": 2500, + "t_step_save": 100, + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 3, + "mapped_weno": 'T', + "weno_eps": 1.0e-6, + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 1, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "rho_wrt": "T", + "parallel_io": "T", + # MHD + "mhd": "T", + "hyper_cleaning": 'T', + "hyper_cleaning_speed": 2.5, + "hyper_cleaning_tau": 0.004, + # Patch 1 - Analytical for v and B + # gamma = 5/3 + # rho = 25/(36*pi) + # p = 5/(12*pi) + # v = (-sin(2*pi*y), sin(2*pi*x), 0) + # B = (-sin(2*pi*y)/sqrt(4*pi), sin(4*pi*x)/sqrt(4*pi), 0) + "patch_icpp(1)%hcid": 250, + "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%length_y": 1.0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 5.0 / (12 * math.pi), + "patch_icpp(1)%Bx": 0.0, + "patch_icpp(1)%By": 0.0, + "patch_icpp(1)%Bz": 0.0, + "patch_icpp(1)%alpha_rho(1)": 25.0 / (36.0 * math.pi), + "patch_icpp(1)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/case_new/2D_orszag_tang_weno5m/case.py b/case_new/2D_orszag_tang_weno5m/case.py new file mode 100644 index 0000000000..29b97ee2d4 --- /dev/null +++ b/case_new/2D_orszag_tang_weno5m/case.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +import json +import math + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0, + "x_domain%end": 1.0, + "y_domain%beg": 0, + "y_domain%end": 1.0, + "m": 512, + "n": 512, + "p": 0, + "dt": 0.0004, + "t_step_start": 0, + "t_step_stop": 2500, + "t_step_save": 100, + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "mapped_weno": 'T', + "weno_eps": 1.0e-6, + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 1, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "rho_wrt": "T", + "parallel_io": "T", + # MHD + "mhd": "T", + "hyper_cleaning": 'T', + "hyper_cleaning_speed": 2.5, + "hyper_cleaning_tau": 0.004, + # Patch 1 - Analytical for v and B + # gamma = 5/3 + # rho = 25/(36*pi) + # p = 5/(12*pi) + # v = (-sin(2*pi*y), sin(2*pi*x), 0) + # B = (-sin(2*pi*y)/sqrt(4*pi), sin(4*pi*x)/sqrt(4*pi), 0) + "patch_icpp(1)%hcid": 250, + "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%length_y": 1.0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 5.0 / (12 * math.pi), + "patch_icpp(1)%Bx": 0.0, + "patch_icpp(1)%By": 0.0, + "patch_icpp(1)%Bz": 0.0, + "patch_icpp(1)%alpha_rho(1)": 25.0 / (36.0 * math.pi), + "patch_icpp(1)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/case_new/2D_orszag_tang_weno5m_nohyper/case.py b/case_new/2D_orszag_tang_weno5m_nohyper/case.py new file mode 100644 index 0000000000..be434674dd --- /dev/null +++ b/case_new/2D_orszag_tang_weno5m_nohyper/case.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +import json +import math + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0, + "x_domain%end": 1.0, + "y_domain%beg": 0, + "y_domain%end": 1.0, + "m": 512, + "n": 512, + "p": 0, + "dt": 0.0004, + "t_step_start": 0, + "t_step_stop": 2500, + "t_step_save": 25, + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "mapped_weno": 'T', + "weno_eps": 1.0e-6, + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 1, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "rho_wrt": "T", + "parallel_io": "T", + # MHD + "mhd": "T", + # "hyper_cleaning": 'T', + # "hyper_cleaning_speed": 2.5, + # "hyper_cleaning_tau": 0.004, + # Patch 1 - Analytical for v and B + # gamma = 5/3 + # rho = 25/(36*pi) + # p = 5/(12*pi) + # v = (-sin(2*pi*y), sin(2*pi*x), 0) + # B = (-sin(2*pi*y)/sqrt(4*pi), sin(4*pi*x)/sqrt(4*pi), 0) + "patch_icpp(1)%hcid": 250, + "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%length_y": 1.0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 5.0 / (12 * math.pi), + "patch_icpp(1)%Bx": 0.0, + "patch_icpp(1)%By": 0.0, + "patch_icpp(1)%Bz": 0.0, + "patch_icpp(1)%alpha_rho(1)": 25.0 / (36.0 * math.pi), + "patch_icpp(1)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/case_new/2D_orszag_tang_weno5z/case.py b/case_new/2D_orszag_tang_weno5z/case.py new file mode 100644 index 0000000000..b51c60956f --- /dev/null +++ b/case_new/2D_orszag_tang_weno5z/case.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +import json +import math + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0, + "x_domain%end": 1.0, + "y_domain%beg": 0, + "y_domain%end": 1.0, + "m": 512, + "n": 512, + "p": 0, + "dt": 0.0004, + "t_step_start": 0, + "t_step_stop": 2500, + "t_step_save": 100, + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "wenoz": 'T', + "weno_eps": 1.0e-6, + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 1, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "rho_wrt": "T", + "parallel_io": "T", + # MHD + "mhd": "T", + "hyper_cleaning": 'T', + "hyper_cleaning_speed": 2.5, + "hyper_cleaning_tau": 0.004, + # Patch 1 - Analytical for v and B + # gamma = 5/3 + # rho = 25/(36*pi) + # p = 5/(12*pi) + # v = (-sin(2*pi*y), sin(2*pi*x), 0) + # B = (-sin(2*pi*y)/sqrt(4*pi), sin(4*pi*x)/sqrt(4*pi), 0) + "patch_icpp(1)%hcid": 250, + "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%length_y": 1.0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 5.0 / (12 * math.pi), + "patch_icpp(1)%Bx": 0.0, + "patch_icpp(1)%By": 0.0, + "patch_icpp(1)%Bz": 0.0, + "patch_icpp(1)%alpha_rho(1)": 25.0 / (36.0 * math.pi), + "patch_icpp(1)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/case_new/2D_orszag_tang_weno7m/case.py b/case_new/2D_orszag_tang_weno7m/case.py new file mode 100644 index 0000000000..83cc2c8fd8 --- /dev/null +++ b/case_new/2D_orszag_tang_weno7m/case.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +import json +import math + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0, + "x_domain%end": 1.0, + "y_domain%beg": 0, + "y_domain%end": 1.0, + "m": 512, + "n": 512, + "p": 0, + "dt": 0.0004, + "t_step_start": 0, + "t_step_stop": 2500, + "t_step_save": 100, + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 7, + "mapped_weno": 'T', + "weno_eps": 1.0e-6, + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 1, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "rho_wrt": "T", + "parallel_io": "T", + # MHD + "mhd": "T", + "hyper_cleaning": 'T', + "hyper_cleaning_speed": 2.5, + "hyper_cleaning_tau": 0.004, + # Patch 1 - Analytical for v and B + # gamma = 5/3 + # rho = 25/(36*pi) + # p = 5/(12*pi) + # v = (-sin(2*pi*y), sin(2*pi*x), 0) + # B = (-sin(2*pi*y)/sqrt(4*pi), sin(4*pi*x)/sqrt(4*pi), 0) + "patch_icpp(1)%hcid": 250, + "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%length_y": 1.0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 5.0 / (12 * math.pi), + "patch_icpp(1)%Bx": 0.0, + "patch_icpp(1)%By": 0.0, + "patch_icpp(1)%Bz": 0.0, + "patch_icpp(1)%alpha_rho(1)": 25.0 / (36.0 * math.pi), + "patch_icpp(1)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) diff --git a/case_new/2D_rotor_weno5m/case.py b/case_new/2D_rotor_weno5m/case.py new file mode 100644 index 0000000000..cdfd6db78d --- /dev/null +++ b/case_new/2D_rotor_weno5m/case.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python3 +import json +import math + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0, + "x_domain%end": 1.0, + "y_domain%beg": 0, + "y_domain%end": 1.0, + "m": 512, + "n": 512, + "p": 0, + "dt": 0.0004, + "t_step_start": 0, + "t_step_stop": 375, + "t_step_save": 25, + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "mapped_weno": 'T', + "weno_eps": 1.0e-6, + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 1, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "rho_wrt": "T", + "parallel_io": "T", + # MHD + "mhd": "T", + "hyper_cleaning": 'T', + "hyper_cleaning_speed": 2.5, + "hyper_cleaning_tau": 0.004, + # Patch 1 - 2D MHD Rotor Problem + # gamma = 1.4 + # Ambient medium (r > 0.1): + # rho = 1 + # p = 1 + # v = (0, 0, 0) + # B = (1, 0, 0) + # Rotor (r <= 0.1): + # rho = 10 + # v has angular velocity of 20 + "patch_icpp(1)%hcid": 252, + "patch_icpp(1)%geometry": 7, + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%length_y": 1.0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%Bx": 1.0, + "patch_icpp(1)%By": 0.0, + "patch_icpp(1)%Bz": 0.0, + "patch_icpp(1)%alpha_rho(1)": 1.0, + "patch_icpp(1)%alpha(1)": 1.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (1.4e00 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + } + ) +) \ No newline at end of file diff --git a/case_new/gauss_div_pulse_2/case.py b/case_new/gauss_div_pulse_2/case.py new file mode 100644 index 0000000000..67a83a0ff3 --- /dev/null +++ b/case_new/gauss_div_pulse_2/case.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +import json + +# ----------------------------------------------------------------------------- +# 2.1 Gaussian Divergence Pulse (1D via uniform y-slice) +# ----------------------------------------------------------------------------- +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain: [0,1]×[0,1] but n=1 → effectively 1D in x + "x_domain%beg":-0.5, + "x_domain%end": 0.5, + "y_domain%beg":-0.5, + "y_domain%end": 0.5, + "m": 100, + "n": 100, + "p": 0, + "dt": 0.005, + "t_step_start": 0, + "t_step_stop": 50, + "t_step_save": 1, + # Numerical Method + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 1, + "weno_order": 1, + "weno_eps": 1.0e-6, + # "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 1, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -1, + "bc_x%end": -1, + "bc_y%beg": -1, + "bc_y%end": -1, + # I/O + # "format": 1, + "format": 2, + "precision": 2, + "prim_vars_wrt": "T", + "rho_wrt": "T", + "parallel_io": "T", + # Physics + "mhd": "T", + "hyper_cleaning": "T", + "hyper_cleaning_speed": 1.5, + "hyper_cleaning_tau": 0.04, + # --- Patch 1: Gaussian Divergence Pulse IC --- + "patch_icpp(1)%geometry": 3, # Cartesian + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%length_x": 1.0, + "patch_icpp(1)%length_y": 1.0, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": 1.0, + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(1)%Bx": 0.0, + "patch_icpp(1)%By": 0.0, + "patch_icpp(1)%Bz": 0.0, + # Fluid EOS + "fluid_pp(1)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + }, + ) +) diff --git a/case_new/slanted_shock_tube/case.py b/case_new/slanted_shock_tube/case.py new file mode 100644 index 0000000000..001bd66d25 --- /dev/null +++ b/case_new/slanted_shock_tube/case.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +import json +import math + +# ----------------------------------------------------------------------------- +# 2.1 Gaussian Divergence Pulse (1D via uniform y-slice) +# ----------------------------------------------------------------------------- +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain: + "x_domain%beg": -1.0, + "x_domain%end": 2.0, + "y_domain%beg":-1.0, + # "y_domain%end": 2.0/256.0, + "y_domain%end": 2.0, + "m": 256*3, + "n": 256*3, + "dt": 0.08/math.sqrt(5.0)/170.0/2.0, + "t_step_start": 0, + "t_step_stop": 170*2, + "t_step_save": 10*2, + # Numerical Method + "num_patches": 1, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 3, + "weno_eps": 1.0e-6, + # "wenoz": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 1, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + "bc_y%beg": -3, + "bc_y%end": -3, + # I/O + # "format": 1, + "format": 2, + "precision": 2, + "prim_vars_wrt": "T", + "rho_wrt": "T", + "parallel_io": "T", + # Physics + "mhd": "T", + # "hyper_cleaning": "T", + # "hyper_cleaning_speed": 30.0, + # "hyper_cleaning_tau": 0.001, + # --- Patch 1: Gaussian Divergence Pulse IC --- + "patch_icpp(1)%hcid": 262, + "patch_icpp(1)%geometry": 7, # Cartesian + "patch_icpp(1)%x_centroid": 0.5, + "patch_icpp(1)%y_centroid": 0.5, + "patch_icpp(1)%length_x": 4.0, + "patch_icpp(1)%length_y": 4.0, + "patch_icpp(1)%vel(1)": 1.0, + "patch_icpp(1)%vel(2)": 1.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 6.0, + "patch_icpp(1)%alpha_rho(1)": 1.0, + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(1)%Bx": 0.0, + "patch_icpp(1)%By": 0.0, + "patch_icpp(1)%Bz": 1/math.sqrt(4*math.pi), + # "patch_icpp(1)%Bz": 1/(4*math.pi), + # "patch_icpp(1)%Bz": 0.0, + # Fluid EOS + "fluid_pp(1)%gamma": 1.0e00 / (5.0 / 3.0 - 1.0e00), + "fluid_pp(1)%pi_inf": 0.0, + }, + ) +) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index dbbf4457e1..8c44d31095 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -1152,6 +1152,8 @@ contains if (cont_damage) qK_prim_vf(damage_idx)%sf(j, k, l) = qK_cons_vf(damage_idx)%sf(j, k, l) + if (hyper_cleaning) qK_prim_vf(psi_idx)%sf(j, k, l) = qK_cons_vf(psi_idx)%sf(j, k, l) + end do end do end do @@ -1420,6 +1422,8 @@ contains if (cont_damage) q_cons_vf(damage_idx)%sf(j, k, l) = q_prim_vf(damage_idx)%sf(j, k, l) + if (hyper_cleaning) q_cons_vf(psi_idx)%sf(j, k, l) = q_prim_vf(psi_idx)%sf(j, k, l) + end do end do end do diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 202c287c34..aa736ff887 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -349,6 +349,9 @@ contains ! Damage state variable if (cont_damage) dbvars = dbvars + 1 + ! Hyperbolic cleaning for MHD + if (hyper_cleaning) dbvars = dbvars + 1 + ! Magnetic field if (mhd) then if (n == 0) then diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 00c5f0ec27..2e657c1824 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -113,6 +113,7 @@ module m_global_parameters integer :: b_size !< Number of components in the b tensor integer :: tensor_size !< Number of components in the nonsymmetric tensor logical :: cont_damage !< Continuum damage modeling + logical :: hyper_cleaning !< Hyperbolic cleaning for MHD logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling !> @} @@ -136,6 +137,7 @@ module m_global_parameters integer :: c_idx !< Index of color function type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns. integer :: damage_idx !< Index of damage state variable (D) for continuum damage model + integer :: psi_idx !< Index of hyperbolic cleaning state variable for MHD !> @} ! Cell Indices for the (local) interior points (O-m, O-n, 0-p). @@ -366,6 +368,7 @@ contains b_size = dflt_int tensor_size = dflt_int cont_damage = .false. + hyper_cleaning = .false. bc_x%beg = dflt_int; bc_x%end = dflt_int bc_y%beg = dflt_int; bc_y%end = dflt_int @@ -725,6 +728,13 @@ contains damage_idx = dflt_int end if + if (hyper_cleaning) then + psi_idx = sys_size + 1 + sys_size = psi_idx + else + psi_idx = dflt_int + end if + end if if (chemistry) then diff --git a/src/post_process/m_mpi_proxy.fpp b/src/post_process/m_mpi_proxy.fpp index 5154d19102..8b67abb312 100644 --- a/src/post_process/m_mpi_proxy.fpp +++ b/src/post_process/m_mpi_proxy.fpp @@ -171,7 +171,7 @@ contains & 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', & & 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', & & 'rkck_adap_dt', 'output_partial_domain', 'relativity', & - & 'cont_damage' ] + & 'cont_damage', 'hyper_cleaning' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 8fece2820b..cf005b6703 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -84,7 +84,7 @@ subroutine s_read_input_file relax_model, cf_wrt, sigma, adv_n, ib, num_ibs, & cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, & cfl_target, surface_tension, bubbles_lagrange, rkck_adap_dt, & - sim_data, hyperelasticity, Bx0, relativity, cont_damage + sim_data, hyperelasticity, Bx0, relativity, cont_damage, hyper_cleaning ! Inquiring the status of the post_process.inp file file_loc = 'post_process.inp' @@ -426,6 +426,14 @@ subroutine s_save_data(t_step, varname, pres, c, H) varname(:) = ' ' end if + if (hyper_cleaning) then + q_sf = q_cons_vf(psi_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end) + write (varname, '(A)') 'psi' + call s_write_variable_to_formatted_database_file(varname, t_step) + + varname(:) = ' ' + end if + ! Adding the pressure to the formatted database file if (pres_wrt .or. prim_vars_wrt) then q_sf = q_prim_vf(E_idx)%sf(x_beg:x_end, y_beg:y_end, z_beg:z_end) diff --git a/src/pre_process/include/2dHardcodedIC.fpp b/src/pre_process/include/2dHardcodedIC.fpp index e4bf547999..18a23e65c7 100644 --- a/src/pre_process/include/2dHardcodedIC.fpp +++ b/src/pre_process/include/2dHardcodedIC.fpp @@ -4,6 +4,19 @@ real(wp) :: r, rmax, gam, umax, p0 real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, intL, alph real(wp) :: factor + real(wp) :: eps_mhd, sigma, C_mhd + real(wp) :: r0, alpha, r2 + real(wp) :: sinA, cosA + + real(wp), parameter :: r_0 = 0.1_wp ! Rotor radius + real(wp), parameter :: r0_sq = r_0*r_0 + real(wp), parameter :: rho_rotor = 10._wp ! Rotor density + real(wp), parameter :: omega = 20._wp ! Angular velocity (v_0/r_0) + real(wp) :: r_sq + + ! Center of the domain + real(wp), parameter :: x_center = 0.5_wp + real(wp), parameter :: y_center = 0.5_wp eps = 1e-9_wp @@ -158,6 +171,91 @@ q_prim_vf(E_idx)%sf(i, j, 0) = 3.e-5_wp end if + ! case 252 is for the 2D MHD Rotor problem + case (252) ! 2D MHD Rotor Problem + ! Ambient conditions are set in the JSON file. + ! This case imposes the dense, rotating cylinder. + ! + ! gamma = 1.4 + ! Ambient medium (r > 0.1): + ! rho = 1, p = 1, v = 0, B = (1,0,0) + ! Rotor (r <= 0.1): + ! rho = 10, p = 1 + ! v has angular velocity w=20, giving v_tan=2 at r=0.1 + + ! Calculate distance squared from the center + r_sq = (x_cc(i) - x_center)**2 + (y_cc(j) - y_center)**2 + + if (r_sq <= r0_sq) then + ! -- Inside the rotor -- + ! Set density + q_prim_vf(contxb)%sf(i, j, 0) = rho_rotor + + ! Set velocity components for rotation + ! v_x = -omega * (y - y_c) + ! v_y = omega * (x - x_c) + q_prim_vf(momxb)%sf(i, j, 0) = -omega * (y_cc(j) - y_center) + q_prim_vf(momxb + 1)%sf(i, j, 0) = omega * (x_cc(i) - x_center) + end if + + case (260) ! Gaussian Divergence Pulse + !--------------------------------------------------------- + ! Bx(x) = 1 + C * erf((x-0.5)/σ) + ! ⇒ ∂Bx/∂x = C * (2/√π) * exp[-((x-0.5)/σ)**2] * (1/σ) + ! Choose C = ε * σ * √π / 2 ⇒ ∂Bx/∂x = ε * exp[-((x-0.5)/σ)**2] + ! ψ is initialized to zero everywhere. + !--------------------------------------------------------- + + eps_mhd = patch_icpp(patch_id)%a(2) + sigma = patch_icpp(patch_id)%a(3) + C_mhd = eps_mhd * sigma * sqrt(pi) * 0.5_wp + + ! B-field + q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp + C_mhd * erf( (x_cc(i) - 0.5_wp) / sigma ) + + case (261) ! Blob + r0 = 1._wp/sqrt(8._wp) + r2 = x_cc(i)**2 + y_cc(j)**2 + r = sqrt(r2) + alpha = r/r0 + if (alpha < 1) then + q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/sqrt(4._wp*pi) * (alpha**8 - 2._wp*alpha**4 + 1._wp) + ! q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/sqrt(4000._wp*pi) * (4096._wp*r2**4 - 128._wp*r2**2 + 1._wp) + ! q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/(4._wp*pi) * (alpha**8 - 2._wp*alpha**4 + 1._wp) + ! q_prim_vf(E_idx)%sf(i,j,0) = 6._wp - q_prim_vf(B_idx%beg)%sf(i,j,0)**2/2._wp + end if + + case (262) ! Tilted 2D MHD shock‐tube at α = arctan2 (≈63.4°) + ! rotate by α = atan(2) + alpha = atan(2._wp) + cosA = cos(alpha) + sinA = sin(alpha) + ! projection along shock normal + r = x_cc(i)*cosA + y_cc(j)*sinA + + if (r <= 0.5_wp) then + ! LEFT state: ρ=1, v∥=+10, v⊥=0, p=20, B∥=B⊥=5/√(4π) + q_prim_vf(contxb)%sf(i,j,0) = 1._wp + q_prim_vf(momxb)%sf(i,j,0) = 10._wp * cosA + q_prim_vf(momxb+1)%sf(i,j,0) = 10._wp * sinA + q_prim_vf(E_idx)%sf(i,j,0) = 20._wp + q_prim_vf(B_idx%beg)%sf(i,j,0) = (5._wp/sqrt(4._wp*pi)) * cosA & + - (5._wp/sqrt(4._wp*pi)) * sinA + q_prim_vf(B_idx%beg+1)%sf(i,j,0) = (5._wp/sqrt(4._wp*pi)) * sinA & + + (5._wp/sqrt(4._wp*pi)) * cosA + else + ! RIGHT state: ρ=1, v∥=−10, v⊥=0, p=1, B∥=B⊥=5/√(4π) + q_prim_vf(contxb)%sf(i,j,0) = 1._wp + q_prim_vf(momxb)%sf(i,j,0) = -10._wp * cosA + q_prim_vf(momxb+1)%sf(i,j,0) = -10._wp * sinA + q_prim_vf(E_idx)%sf(i,j,0) = 1._wp + q_prim_vf(B_idx%beg)%sf(i,j,0) = (5._wp/sqrt(4._wp*pi)) * cosA & + - (5._wp/sqrt(4._wp*pi)) * sinA + q_prim_vf(B_idx%beg+1)%sf(i,j,0) = (5._wp/sqrt(4._wp*pi)) * sinA & + + (5._wp/sqrt(4._wp*pi)) * cosA + end if + ! v^z and B^z remain zero by default + case default if (proc_rank == 0) then call s_int_to_str(patch_id, iStr) diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index cec3227e2e..8c817e5560 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -94,6 +94,7 @@ module m_global_parameters integer :: tensor_size !< Number of components in the nonsymmetric tensor logical :: pre_stress !< activate pre_stressed domain logical :: cont_damage !< continuum damage modeling + logical :: hyper_cleaning !< Hyperbolic cleaning for MHD logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling ! Annotations of the structure, i.e. the organization, of the state vectors @@ -113,6 +114,7 @@ module m_global_parameters integer :: c_idx !< Index of the color function type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns. integer :: damage_idx !< Index of damage state variable (D) for continuum damage model + integer :: psi_idx !< Index of hyperbolic cleaning state variable for MHD ! Cell Indices for the (local) interior points (O-m, O-n, 0-p). ! Stands for "InDices With BUFFer". @@ -343,6 +345,7 @@ contains b_size = dflt_int tensor_size = dflt_int cont_damage = .false. + hyper_cleaning = .false. mhd = .false. relativity = .false. @@ -816,6 +819,11 @@ contains sys_size = damage_idx end if + if (hyper_cleaning) then + psi_idx = sys_size + 1 + sys_size = psi_idx + end if + end if if (chemistry) then diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index dd2083b805..47c7f96d28 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -119,6 +119,13 @@ contains q_prim_vf(damage_idx)%sf = 0._wp end if + ! Initial hyper_cleaning state is always zero + ! TODO more general + if (hyper_cleaning) then + q_cons_vf(psi_idx)%sf = 0._wp + q_prim_vf(psi_idx)%sf = 0._wp + end if + ! Setting default values for patch identities bookkeeping variable. ! This is necessary to avoid any confusion in the assessment of the ! extent of application that the overwrite permissions give a patch diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index c9888e79ca..51019c8c19 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -57,7 +57,8 @@ contains & 'qbmm', 'file_per_process', 'adv_n', 'ib' , 'cfl_adap_dt', & & 'cfl_const_dt', 'cfl_dt', 'surface_tension', & & 'hyperelasticity', 'pre_stress', 'elliptic_smoothing', 'viscous',& - & 'bubbles_lagrange', 'mhd', 'relativity', 'cont_damage' ] + & 'bubbles_lagrange', 'mhd', 'relativity', 'cont_damage', & + & 'hyper_cleaning' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index ee2ed3f30e..e36d5cad0b 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -144,7 +144,7 @@ contains sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, & n_start_old, surface_tension, hyperelasticity, pre_stress, & rkck_adap_dt, elliptic_smoothing, elliptic_smoothing_iters, & - viscous, bubbles_lagrange, Bx0, relativity, cont_damage + viscous, bubbles_lagrange, Bx0, relativity, cont_damage, hyper_cleaning ! Inquiring the status of the pre_process.inp file file_loc = 'pre_process.inp' @@ -822,6 +822,8 @@ contains real(wp), intent(inout) :: time_avg, time_final logical, intent(inout) :: file_exists + integer :: j, k + ! Setting up the grid and the initial condition. If the grid is read in from ! preexisting grid data files, it is checked for consistency. If the grid is ! not read in, it is generated from scratch according to the inputs provided @@ -838,6 +840,16 @@ contains call s_generate_initial_condition() + ! hard-coded psi + if (hyper_cleaning) then + do j = 0, m + do k = 0, n + q_cons_vf(psi_idx)%sf(j,k,0) = 1d-2*exp(-( x_cc(j)**2 + y_cc(k)**2 )/(2.0*0.05**2)) + q_prim_vf(psi_idx)%sf(j,k,0) = q_cons_vf(psi_idx)%sf(j,k,0) + end do + end do + end if + if (relax) then if (proc_rank == 0) then print *, 'initial condition might have been altered due to enforcement of & diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 56b21207ff..9e3045f95f 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1440,6 +1440,11 @@ contains tmp = damage_state call s_mpi_allreduce_sum(tmp, damage_state) end if + + ! if (hyper_cleaning) then + ! tmp = 0._wp + ! call s_mpi_allreduce_sum(tmp, tmp) ! TODO + ! end if end if if (proc_rank == 0) then if (n == 0) then diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 329153bfdf..4ab6259c50 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -157,6 +157,7 @@ module m_global_parameters logical :: shear_stress !< Shear stresses logical :: bulk_stress !< Bulk stresses logical :: cont_damage !< Continuum damage modeling + logical :: hyper_cleaning !< Hyperbolic cleaning for MHD for divB=0 !$acc declare create(chemistry) @@ -177,7 +178,7 @@ module m_global_parameters !$acc declare create(num_dims, num_vels, weno_polyn, weno_order, weno_num_stencils, num_fluids, wenojs, mapped_weno, wenoz, teno, wenoz_q, mhd, relativity) #:endif - !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, hyperelasticity, hyper_model, elasticity, low_Mach, viscous, shear_stress, bulk_stress, cont_damage) + !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, hyperelasticity, hyper_model, elasticity, low_Mach, viscous, shear_stress, bulk_stress, cont_damage, hyper_cleaning) logical :: relax !< activate phase change integer :: relax_model !< Relaxation model @@ -243,6 +244,7 @@ module m_global_parameters type(int_bounds_info) :: species_idx !< Indexes of first & last concentration eqns. integer :: c_idx !< Index of color function integer :: damage_idx !< Index of damage state variable (D) for continuum damage model + integer :: psi_idx !< Index of hyperbolic cleaning state variable for MHD !> @} !$acc declare create(bub_idx) @@ -497,6 +499,13 @@ module m_global_parameters !$acc declare create(tau_star, cont_damage_s, alpha_bar) !> @} + !> @name MHD Hyperbolic cleaning parameters + !> @{! + real(wp) :: hyper_cleaning_speed !< Hyperbolic cleaning wave speed (c_h) + real(wp) :: hyper_cleaning_tau !< Hyperbolic cleaning tau + !$acc declare create(hyper_cleaning_speed, hyper_cleaning_tau) + !> @} + contains !> Assigns default values to the user inputs before reading @@ -568,6 +577,7 @@ contains shear_stress = .false. bulk_stress = .false. cont_damage = .false. + hyper_cleaning = .false. #:if not MFC_CASE_OPTIMIZATION mapped_weno = .false. @@ -761,6 +771,8 @@ contains ! MHD Bx0 = dflt_real powell = .false. + hyper_cleaning_speed = dflt_real + hyper_cleaning_tau = dflt_real #:if not MFC_CASE_OPTIMIZATION mhd = .false. @@ -1104,6 +1116,11 @@ contains sys_size = damage_idx end if + if (hyper_cleaning) then + psi_idx = sys_size + 1 + sys_size = psi_idx + end if + end if ! END: Volume Fraction Model @@ -1225,6 +1242,8 @@ contains !$acc update device(cont_damage, tau_star, cont_damage_s, alpha_bar) + !$acc update device(hyper_cleaning, hyper_cleaning_speed, hyper_cleaning_tau) + #:if not MFC_CASE_OPTIMIZATION !$acc update device(wenojs, mapped_weno, wenoz, teno) !$acc update device(wenoz_q) diff --git a/src/simulation/m_mhd.fpp b/src/simulation/m_mhd.fpp index 2a188fc332..bba4e73460 100644 --- a/src/simulation/m_mhd.fpp +++ b/src/simulation/m_mhd.fpp @@ -17,8 +17,8 @@ module m_mhd implicit none - private; public :: s_initialize_mhd_powell_module, & - s_finalize_mhd_powell_module, & + private; public :: s_initialize_mhd_clean_module, & + s_finalize_mhd_clean_module, & s_compute_mhd_powell_rhs real(wp), allocatable, dimension(:, :, :) :: du_dx, du_dy, du_dz @@ -33,7 +33,7 @@ module m_mhd contains - subroutine s_initialize_mhd_powell_module + subroutine s_initialize_mhd_clean_module integer :: i, k, r @@ -62,17 +62,19 @@ contains !$acc update device(fd_coeff_z_h) end if - end subroutine s_initialize_mhd_powell_module + end subroutine s_initialize_mhd_clean_module !> Compute the Powell source term to correct the magnetic field divergence. !! The Powell source term is: !! S = - (divB) [ 0, Bx, By, Bz, vdotB, vx, vy, vz ]^T !! @param q_prim_vf Primitive variables !! @param rhs_vf rhs variables - subroutine s_compute_mhd_powell_rhs(q_prim_vf, rhs_vf) + subroutine s_compute_mhd_powell_rhs(idir, q_prim_vf, rhs_vf, flux_gsrc_n) + integer, intent(in) :: idir type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf + type(scalar_field), dimension(sys_size), intent(in) :: flux_gsrc_n integer :: k, l, q, r real(wp), dimension(3) :: v, B @@ -84,20 +86,15 @@ contains do l = 0, n do k = 0, m - divB = 0._wp - !$acc loop seq - do r = -fd_number, fd_number - divB = divB + q_prim_vf(B_idx%beg)%sf(k + r, l, q)*fd_coeff_x_h(r, k) - end do - !$acc loop seq - do r = -fd_number, fd_number - divB = divB + q_prim_vf(B_idx%beg + 1)%sf(k, l + r, q)*fd_coeff_y_h(r, l) - end do - if (p > 0) then - !$acc loop seq - do r = -fd_number, fd_number - divB = divB + q_prim_vf(B_idx%beg + 2)%sf(k, l, q + r)*fd_coeff_z_h(r, q) - end do + if (idir == 1) then + divB = 1._wp/dx(k)* & + (flux_gsrc_n(1)%sf(k, l, q) - flux_gsrc_n(1)%sf(k - 1, l, q)) + else if (idir == 2) then + divB = 1._wp/dy(l)* & + (flux_gsrc_n(1)%sf(k, l, q) - flux_gsrc_n(1)%sf(k, l - 1, q)) + else if (idir == 3) then + divB = 1._wp/dz(l)* & + (flux_gsrc_n(1)%sf(k, l, q) - flux_gsrc_n(1)%sf(k, l, q - 1)) end if v(1) = q_prim_vf(momxb)%sf(k, l, q) @@ -136,7 +133,7 @@ contains end subroutine s_compute_mhd_powell_rhs - subroutine s_finalize_mhd_powell_module + subroutine s_finalize_mhd_clean_module @:DEALLOCATE(du_dx, dv_dx, dw_dx) @:DEALLOCATE(fd_coeff_x_h) @@ -147,6 +144,6 @@ contains @:DEALLOCATE(fd_coeff_z_h) end if - end subroutine s_finalize_mhd_powell_module + end subroutine s_finalize_mhd_clean_module end module m_mhd diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 9dd922f36b..596dbd94ed 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -136,7 +136,8 @@ contains & 'bc_z%grcbc_in', 'bc_z%grcbc_out', 'bc_z%grcbc_vel_out', & & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'surface_tension', & & 'viscous', 'shear_stress', 'bulk_stress', 'bubbles_lagrange', & - & 'hyperelasticity', 'rkck_adap_dt', 'powell', 'cont_damage' ] + & 'hyperelasticity', 'rkck_adap_dt', 'powell', 'cont_damage', & + & 'hyper_cleaning' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -175,7 +176,8 @@ contains & 'x_domain%beg', 'x_domain%end', 'y_domain%beg', 'y_domain%end', & & 'z_domain%beg', 'z_domain%end', 'x_a', 'x_b', 'y_a', 'y_b', 'z_a', & & 'z_b', 't_stop', 't_save', 'cfl_target', 'rkck_tolerance', 'Bx0', & - & 'tau_star', 'cont_damage_s', 'alpha_bar' ] + & 'tau_star', 'cont_damage_s', 'alpha_bar', 'hyper_cleaning_speed', & + & 'hyper_cleaning_tau' ] call MPI_BCAST(${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #:endfor diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 09e4a17390..0fcfca3c41 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -190,7 +190,7 @@ contains @:ALLOCATE(q_prim_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end)) end do - num_eqns_after_adv = count((/surface_tension, cont_damage/)) + num_eqns_after_adv = count((/surface_tension, cont_damage, hyper_cleaning/)) do l = adv_idx%end + 1, sys_size - num_eqns_after_adv @:ALLOCATE(q_prim_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end)) @@ -229,6 +229,13 @@ contains !$acc enter data attach(q_prim_qp%vf(damage_idx)%sf) end if + if (hyper_cleaning) then + q_prim_qp%vf(psi_idx)%sf => & + q_cons_qp%vf(psi_idx)%sf + !$acc enter data copyin(q_prim_qp%vf(psi_idx)%sf) + !$acc enter data attach(q_prim_qp%vf(psi_idx)%sf) + end if + if (viscous) then @:ALLOCATE(tau_Re_vf(1:sys_size)) do i = 1, num_dims @@ -856,7 +863,7 @@ contains ! END: Additional physics and source terms call nvtxStartRange("RHS-MHD") - if (mhd .and. powell) call s_compute_mhd_powell_rhs(q_prim_qp%vf, rhs_vf) + if (powell) call s_compute_mhd_powell_rhs(id, q_prim_qp%vf, rhs_vf, flux_gsrc_n(id)%vf) call nvtxEndRange end do diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 69a8b4a61b..9effe5f18e 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -677,6 +677,11 @@ contains call s_compute_fast_magnetosonic_speed(rho_R, c_R, B%R, norm_dir, c_fast%R, H_R) end if + if (hyper_cleaning) then ! mhd + c_fast%L = min(c_fast%L, -hyper_cleaning_speed) + c_fast%R = max(c_fast%R, hyper_cleaning_speed) + end if + if (viscous) then !$acc loop seq do i = 1, 2 @@ -1035,28 +1040,39 @@ contains else ! 2D/3D: Bx, By, Bz /= const. but zero flux component in the same direction ! B_x d/d${XYZ}$ flux = (1 - delta(x,${XYZ}$)) * (v_${XYZ}$ * B_x - v_x * B_${XYZ}$) - flux_rs${XYZ}$_vf(j, k, l, B_idx%beg) = (1 - dir_flg(1))*( & + flux_rs${XYZ}$_vf(j, k, l, B_idx%beg) = ( & s_M*(vel_R(dir_idx(1))*B%R(1) - vel_R(1)*B%R(norm_dir)) - & s_P*(vel_L(dir_idx(1))*B%L(1) - vel_L(1)*B%L(norm_dir)) + & s_M*s_P*(B%L(1) - B%R(1)))/(s_M - s_P) ! B_y d/d${XYZ}$ flux = (1 - delta(y,${XYZ}$)) * (v_${XYZ}$ * B_y - v_y * B_${XYZ}$) - flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1) = (1 - dir_flg(2))*( & + flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1) = ( & s_M*(vel_R(dir_idx(1))*B%R(2) - vel_R(2)*B%R(norm_dir)) - & s_P*(vel_L(dir_idx(1))*B%L(2) - vel_L(2)*B%L(norm_dir)) + & s_M*s_P*(B%L(2) - B%R(2)))/(s_M - s_P) ! B_z d/d${XYZ}$ flux = (1 - delta(z,${XYZ}$)) * (v_${XYZ}$ * B_z - v_z * B_${XYZ}$) - flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + 2) = (1 - dir_flg(3))*( & + flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + 2) = ( & s_M*(vel_R(dir_idx(1))*B%R(3) - vel_R(3)*B%R(norm_dir)) - & s_P*(vel_L(dir_idx(1))*B%L(3) - vel_L(3)*B%L(norm_dir)) + & s_M*s_P*(B%L(3) - B%R(3)))/(s_M - s_P) + if (hyper_cleaning) then + flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + norm_dir - 1) = flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + norm_dir - 1) + & + (s_M*qR_prim_rs${XYZ}$_vf(j + 1, k, l, psi_idx) - s_P*qL_prim_rs${XYZ}$_vf(j, k, l, psi_idx))/(s_M - s_P) + + flux_rs${XYZ}$_vf(j, k, l, psi_idx) = (hyper_cleaning_speed**2*(s_M*B%R(norm_dir) - s_P*B%L(norm_dir)) + s_M*s_P*(qL_prim_rs${XYZ}$_vf(j, k, l, psi_idx) - qR_prim_rs${XYZ}$_vf(j + 1, k, l, psi_idx)))/(s_M - s_P) + end if + end if flux_src_rs${XYZ}$_vf(j, k, l, advxb) = 0._wp end if + if (powell) then + flux_gsrc_rs${XYZ}$_vf(j, k, l, 1) = flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + norm_dir - 1) + end if + #:if (NORM_DIR == 2) if (cyl_coord) then !Substituting the advective flux into the inviscid geometrical source flux @@ -4951,7 +4967,7 @@ contains end do end do - if (cyl_coord) then + if (cyl_coord .or. powell) then !$acc parallel loop collapse(4) gang vector default(present) do i = 1, sys_size do l = is3%beg, is3%end @@ -5003,7 +5019,7 @@ contains end do end do end do - if (grid_geometry == 3) then + if (grid_geometry == 3 .or. powell) then !$acc parallel loop collapse(4) gang vector default(present) do i = 1, sys_size do j = is1%beg, is1%end @@ -5055,6 +5071,20 @@ contains end do end do + if (powell) then + !$acc parallel loop collapse(4) gang vector default(present) + do i = 1, sys_size + do l = is3%beg, is3%end + do k = is2%beg, is2%end + do j = is1%beg, is1%end + flux_gsrc_vf(i)%sf(j, k, l) = & + flux_gsrc_rsx_vf(j, k, l, i) + end do + end do + end do + end do + end if + !$acc parallel loop collapse(3) gang vector default(present) do l = is3%beg, is3%end do k = is2%beg, is2%end diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 0606ad5ce1..d4e64f57d6 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -176,7 +176,8 @@ contains bubbles_lagrange, lag_params, & rkck_adap_dt, rkck_tolerance, & hyperelasticity, R0ref, Bx0, powell, & - cont_damage, tau_star, cont_damage_s, alpha_bar + cont_damage, tau_star, cont_damage_s, alpha_bar, & + hyper_cleaning, hyper_cleaning_speed, hyper_cleaning_tau ! Checking that an input file has been provided by the user. If it ! has, then the input file is read in, otherwise, simulation exits. @@ -1528,7 +1529,7 @@ contains if (hypoelasticity) call s_initialize_hypoelastic_module() if (hyperelasticity) call s_initialize_hyperelastic_module() - if (mhd .and. powell) call s_initialize_mhd_powell_module + if (powell) call s_initialize_mhd_clean_module end subroutine s_initialize_modules @@ -1667,7 +1668,7 @@ contains if (surface_tension) call s_finalize_surface_tension_module() if (bodyForces) call s_finalize_body_forces_module() - if (mhd .and. powell) call s_finalize_mhd_powell_module + if (powell) call s_finalize_mhd_clean_module ! Terminating MPI execution environment call s_mpi_finalize() diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 045dfda597..0fd84d28d7 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -188,6 +188,13 @@ contains @:ACC_SETUP_SFs(q_prim_vf(damage_idx)) end if + if (hyper_cleaning) then + @:ALLOCATE(q_prim_vf(psi_idx)%sf(idwbuff(1)%beg:idwbuff(1)%end, & + idwbuff(2)%beg:idwbuff(2)%end, & + idwbuff(3)%beg:idwbuff(3)%end)) + @:ACC_SETUP_SFs(q_prim_vf(psi_idx)) + end if + if (model_eqns == 3) then do i = internalEnergies_idx%beg, internalEnergies_idx%end @:ALLOCATE(q_prim_vf(i)%sf(idwbuff(1)%beg:idwbuff(1)%end, & @@ -421,6 +428,21 @@ contains end if end if + if (hyper_cleaning) then + do l = 0, p + do k = 0, n + do j = 0, m + ! q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & + ! q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) & + ! - dt*q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)/hyper_cleaning_tau + q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & + q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)* & + exp(-dt/hyper_cleaning_tau) + end do + end do + end do + end if + call nvtxEndRange end subroutine s_1st_order_tvd_rk @@ -716,6 +738,21 @@ contains end if end if + if (hyper_cleaning) then + do l = 0, p + do k = 0, n + do j = 0, m + ! q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & + ! q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) & + ! - dt*q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)/hyper_cleaning_tau + q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & + q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)* & + exp(-dt/hyper_cleaning_tau) + end do + end do + end do + end if + ! Stage 2 of 3 call s_compute_rhs(q_cons_ts(2)%vf, q_T_sf, q_prim_vf, rhs_vf, pb_ts(2)%sf, rhs_pb, mv_ts(2)%sf, rhs_mv, t_step, time_avg) @@ -793,6 +830,21 @@ contains end if end if + if (hyper_cleaning) then + do l = 0, p + do k = 0, n + do j = 0, m + ! q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & + ! q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) & + ! - dt*q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)/hyper_cleaning_tau + q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & + q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)* & + exp(-dt/hyper_cleaning_tau) + end do + end do + end do + end if + ! Stage 3 of 3 call s_compute_rhs(q_cons_ts(2)%vf, q_T_sf, q_prim_vf, rhs_vf, pb_ts(2)%sf, rhs_pb, mv_ts(2)%sf, rhs_mv, t_step, time_avg) @@ -873,6 +925,21 @@ contains end if end if + if (hyper_cleaning) then + do l = 0, p + do k = 0, n + do j = 0, m + ! q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & + ! q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) & + ! - dt*q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)/hyper_cleaning_tau + q_cons_ts(1)%vf(psi_idx)%sf(j, k, l) = & + q_cons_ts(1)%vf(psi_idx)%sf(j, k, l)* & + exp(-dt/hyper_cleaning_tau) + end do + end do + end do + end if + if (.not. adap_dt) then call nvtxEndRange call cpu_time(finish) @@ -1255,6 +1322,10 @@ contains @:DEALLOCATE(q_prim_vf(damage_idx)%sf) end if + if (hyper_cleaning) then + @:DEALLOCATE(q_prim_vf(psi_idx)%sf) + end if + if (bubbles_euler) then do i = bub_idx%beg, bub_idx%end @:DEALLOCATE(q_prim_vf(i)%sf) diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 0294ffc5e5..10bcf7d087 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -61,6 +61,7 @@ def analytic(self): 'Bx0': ParamType.REAL, 'relativity': ParamType.LOG, 'cont_damage': ParamType.LOG, + 'hyper_cleaning': ParamType.LOG, } PRE_PROCESS = COMMON.copy() @@ -260,6 +261,8 @@ def analytic(self): 'tau_star': ParamType.REAL, 'cont_damage_s': ParamType.REAL, 'alpha_bar': ParamType.REAL, + 'hyper_cleaning_speed': ParamType.REAL, + 'hyper_cleaning_tau': ParamType.REAL, }) for var in [ 'heatTransfer_model', 'massTransfer_model', 'pressure_corrector',