diff --git a/.github/workflows/lint-source.yml b/.github/workflows/lint-source.yml index f9f88228e7..fa2d300943 100644 --- a/.github/workflows/lint-source.yml +++ b/.github/workflows/lint-source.yml @@ -42,6 +42,10 @@ jobs: run: | ! grep -iR -e '\.\.\.' -e '\-\-\-' -e '===' ./src/* + - name: Looking for all caps comments in source + run: | + ! grep -RE '^\s+[!]\s+[A-Z]{5}' ./src/* + - name: Looking for junk comments in examples run: | ! grep -R '# ===' ./benchmarks **/*.py @@ -49,3 +53,4 @@ jobs: ! grep -R '===' ./benchmarks/**/*.py ! grep -R '===' ./examples/**/*.py + diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 07ee6bb282..79a24221f6 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -523,6 +523,7 @@ To restart the simulation from $k$-th time step, see [Restarting Cases](running. | `vel_wrt(i)` | Logical | Add the $i$-direction velocity to the database | | `E_wrt` | Logical | Add the total energy to the database | | `pres_wrt` | Logical | Add the pressure to the database | +| `tau_wrt` | Logical | Add the elastic stresses to the database | | `alpha_wrt(i)` | Logical | Add the volume fraction of fluid $i$ to the database | | `gamma_wrt` | Logical | Add the specific heat ratio function to the database | | `heat_ratio_wrt` | Logical | Add the specific heat ratio to the database | diff --git a/examples/1D_hyper_impact_strong/case.py b/examples/1D_hyper_impact_strong/case.py new file mode 100755 index 0000000000..de9b1effff --- /dev/null +++ b/examples/1D_hyper_impact_strong/case.py @@ -0,0 +1,92 @@ +#!/usr/bin/env python3 +import math +import json + +# Numerical setup +Nx = 201 +dx = 1.0 / (1.0 * (Nx + 1)) + +Tend = 41e-06 +Nt = 4000 +mydt = Tend / (1.0 * Nt) + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0.0e00, + "x_domain%end": 1.0e00, + # 'y_domain%beg' : 0.E+00, + # 'y_domain%end' : 0.002, + "m": Nx, + "n": 0, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": int(Nt), + "t_step_save": int(Nt / 200), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F", + "null_weights": "F", + "mp_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + #'bc_y%beg' : -3, + #'bc_y%end' : -3, + # Turning on Hypoelasticity + "hyperelasticity": "T", + "hyper_model": 1, + "fd_order": 4, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + # Patch 1 L + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.25, + # 'patch_icpp(1)%y_centroid' : 0.001, + "patch_icpp(1)%length_x": 0.5, + # 'patch_icpp(1)%length_y' : 0.002, + "patch_icpp(1)%vel(1)": 1000, + # 'patch_icpp(1)%vel(2)' : 100*0, + "patch_icpp(1)%pres": 1.0e5, + "patch_icpp(1)%alpha_rho(1)": 1000, + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(1)%tau_e(1)": 0.0, + # Patch 2 R + "patch_icpp(2)%geometry": 1, + "patch_icpp(2)%x_centroid": 0.75, + # 'patch_icpp(2)%y_centroid' : 0.001, + "patch_icpp(2)%length_x": 0.5, + # 'patch_icpp(2)%length_y' : 0.002, + "patch_icpp(2)%vel(1)": 1000, + # 'patch_icpp(2)%vel(2)' : -100*0, + "patch_icpp(2)%pres": 1.0e05, + "patch_icpp(2)%alpha_rho(1)": 1000, + "patch_icpp(2)%alpha(1)": 1.0, + "patch_icpp(2)%tau_e(1)": 0.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00), + "fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00), + "fluid_pp(1)%G": 1e010, + } + ) +) diff --git a/examples/1D_hyper_impact_weak/case.py b/examples/1D_hyper_impact_weak/case.py new file mode 100755 index 0000000000..87a24c330b --- /dev/null +++ b/examples/1D_hyper_impact_weak/case.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 +import math +import json + +# Numerical setup +Nx = 201 +dx = 1.0 / (1.0 * (Nx + 1)) + +Tend = 64e-06 +Nt = 4000 +mydt = Tend / (1.0e00 * Nt) +# print(mydt) +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0.0e00, + "x_domain%end": 1.0e00, + "m": Nx, + "n": 0, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": int(Nt), + "t_step_save": int(Nt / 200), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 3, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F", + "null_weights": "F", + "mp_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + # Turning on Hyperelasticity + "hyperelasticity": "T", + "hyper_model": 1, + "fd_order": 4, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + # Patch 1 L + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.25, + "patch_icpp(1)%length_x": 0.5, + "patch_icpp(1)%vel(1)": 10, + "patch_icpp(1)%pres": 1.0e5, + "patch_icpp(1)%alpha_rho(1)": 1000, + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(1)%tau_e(1)": 0.0, + # Patch 2 R + "patch_icpp(2)%geometry": 1, + "patch_icpp(2)%x_centroid": 0.75, + "patch_icpp(2)%length_x": 0.5, + "patch_icpp(2)%vel(1)": -10, # 10, + "patch_icpp(2)%pres": 1.0e05, + "patch_icpp(2)%alpha_rho(1)": 1000, + "patch_icpp(2)%alpha(1)": 1.0, + "patch_icpp(2)%tau_e(1)": 0.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00), + "fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00), + "fluid_pp(1)%G": 1.0e00, # .E+010, + } + ) +) diff --git a/examples/1D_hyper_impact_weak/original.py b/examples/1D_hyper_impact_weak/original.py new file mode 100755 index 0000000000..3ac4a61404 --- /dev/null +++ b/examples/1D_hyper_impact_weak/original.py @@ -0,0 +1,90 @@ +#!/usr/bin/env python3 +import math +import json + +# Numerical setup +Nx = 201 +dx = 1.0 / (1.0 * (Nx + 1)) + +Tend = 64e-06 +Nt = 4000 +mydt = Tend / (1.0e00 * Nt) +# print(mydt) +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0.0e00, + "x_domain%end": 1.0e00, + # 'y_domain%beg' : 0.E+00, + # 'y_domain%end' : 0.002, + "m": Nx, + "n": 0, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": int(Nt), + "t_step_save": int(Nt / 200), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 3, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F", + "null_weights": "F", + "mp_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + #'bc_y%beg' : -3, + #'bc_y%end' : -3, + # Turning on Hyperelasticity + "hyperelasticity": "T", + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + # Patch 1 L + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.25, + # 'patch_icpp(1)%y_centroid' : 0.001, + "patch_icpp(1)%length_x": 0.5, + # 'patch_icpp(1)%length_y' : 0.002, + "patch_icpp(1)%vel(1)": 10, + # 'patch_icpp(1)%vel(2)' : 100*0, + "patch_icpp(1)%pres": 1.0e5, + "patch_icpp(1)%alpha_rho(1)": 1000, + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(1)%tau_e(1)": 0.0, + # Patch 2 R + "patch_icpp(2)%geometry": 1, + "patch_icpp(2)%x_centroid": 0.75, + # 'patch_icpp(2)%y_centroid' : 0.001, + "patch_icpp(2)%length_x": 0.5, + # 'patch_icpp(2)%length_y' : 0.002, + "patch_icpp(2)%vel(1)": 10, # -10, + # 'patch_icpp(2)%vel(2)' : -100*0, + "patch_icpp(2)%pres": 1.0e05, + "patch_icpp(2)%alpha_rho(1)": 1000, + "patch_icpp(2)%alpha(1)": 1.0, + "patch_icpp(2)%tau_e(1)": 0.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00), + "fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00), + "fluid_pp(1)%G": 1.0e00, # .E+010, + } + ) +) diff --git a/examples/1D_hypo_impact_strong/case.py b/examples/1D_hypo_impact_strong/case.py new file mode 100755 index 0000000000..6044fc5cf5 --- /dev/null +++ b/examples/1D_hypo_impact_strong/case.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 +import math +import json + +# Numerical setup +Nx = 201 +dx = 1.0 / (1.0 * (Nx + 1)) + +Tend = 41e-06 +Nt = 4000 +mydt = Tend / (1.0 * Nt) + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0.0e00, + "x_domain%end": 1.0e00, + # 'y_domain%beg' : 0.E+00, + # 'y_domain%end' : 0.002, + "m": Nx, + "n": 0, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": int(Nt), + "t_step_save": int(Nt / 200), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F", + "null_weights": "F", + "mp_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + #'bc_y%beg' : -3, + #'bc_y%end' : -3, + # Turning on Hypoelasticity + "hypoelasticity": "T", + "fd_order": 4, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + # Patch 1 L + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.25, + # 'patch_icpp(1)%y_centroid' : 0.001, + "patch_icpp(1)%length_x": 0.5, + # 'patch_icpp(1)%length_y' : 0.002, + "patch_icpp(1)%vel(1)": 1000, + # 'patch_icpp(1)%vel(2)' : 100*0, + "patch_icpp(1)%pres": 1.0e5, + "patch_icpp(1)%alpha_rho(1)": 1000, + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(1)%tau_e(1)": 0.0, + # Patch 2 R + "patch_icpp(2)%geometry": 1, + "patch_icpp(2)%x_centroid": 0.75, + # 'patch_icpp(2)%y_centroid' : 0.001, + "patch_icpp(2)%length_x": 0.5, + # 'patch_icpp(2)%length_y' : 0.002, + "patch_icpp(2)%vel(1)": 1000, + # 'patch_icpp(2)%vel(2)' : -100*0, + "patch_icpp(2)%pres": 1.0e05, + "patch_icpp(2)%alpha_rho(1)": 1000, + "patch_icpp(2)%alpha(1)": 1.0, + "patch_icpp(2)%tau_e(1)": 0.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00), + "fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00), + "fluid_pp(1)%G": 1e010, + } + ) +) diff --git a/examples/1D_hypo_impact_weak/case.py b/examples/1D_hypo_impact_weak/case.py new file mode 100755 index 0000000000..05d1f5ca94 --- /dev/null +++ b/examples/1D_hypo_impact_weak/case.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +import math +import json + +# Numerical setup +Nx = 201 +dx = 1.0 / (1.0 * (Nx + 1)) + +Tend = 64e-06 +Nt = 4000 +mydt = Tend / (1.0 * Nt) + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": 0.0e00, + "x_domain%end": 1.0e00, + "m": Nx, + "n": 0, + "p": 0, + "dt": mydt, + "t_step_start": 0, + "t_step_stop": int(Nt), + "t_step_save": int(Nt / 200), + # Simulation Algorithm Parameters + "num_patches": 2, + "model_eqns": 2, + "alt_soundspeed": "F", + "num_fluids": 2, + "mpp_lim": "F", + "mixture_err": "F", + "time_stepper": 3, + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "F", + "null_weights": "F", + "mp_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + # Turning on Hypoelasticity + "hypoelasticity": "T", + "fd_order": 4, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + # Patch 1 L + "patch_icpp(1)%geometry": 1, + "patch_icpp(1)%x_centroid": 0.25, + "patch_icpp(1)%length_x": 0.5, + "patch_icpp(1)%vel(1)": 10, + "patch_icpp(1)%pres": 1.0e5, + "patch_icpp(1)%alpha_rho(1)": 1000, + "patch_icpp(1)%alpha_rho(2)": 0.0, + "patch_icpp(1)%alpha(1)": 1.0, + "patch_icpp(1)%alpha(2)": 0.0, + "patch_icpp(1)%tau_e(1)": 0.0, + # Patch 2 R + "patch_icpp(2)%geometry": 1, + "patch_icpp(2)%x_centroid": 0.75, + "patch_icpp(2)%length_x": 0.5, + "patch_icpp(2)%vel(1)": -10, + "patch_icpp(2)%pres": 1.0e05, + "patch_icpp(2)%alpha_rho(1)": 0.0, + "patch_icpp(2)%alpha_rho(2)": 1000, + "patch_icpp(2)%alpha(1)": 0.0, + "patch_icpp(2)%alpha(2)": 1.0, + "patch_icpp(2)%tau_e(1)": 0.0, + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00), + "fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00), + "fluid_pp(1)%G": 1e010, + "fluid_pp(2)%gamma": 1.0e00 / (4.4e00 - 1.0e00), + "fluid_pp(2)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00), + "fluid_pp(2)%G": 1e010, + } + ) +) diff --git a/examples/3D_hyper_bubingel/case.py b/examples/3D_hyper_bubingel/case.py new file mode 100755 index 0000000000..740abc3b93 --- /dev/null +++ b/examples/3D_hyper_bubingel/case.py @@ -0,0 +1,419 @@ +#!/usr/bin/env python3 +import math, json + +## 1 FOR BACKGROUND, 2 FOR BUBBLE, 3 FOR GEL +# Pressure [Pa] +p01 = 5e6 +p02 = 3550 +p03 = p01 + +# Temperature [K] +T01 = 298.15 +T02 = 298.15 +T03 = T01 + +#### FLUID PROPERTIES #### + +### liquid water ### +# pi infty +piwl = 1.0e09 +# qv +qvwl = -1167000 +# qv' +qvpwl = 0.0e0 +# cv +cvwl = 1816 +# cp +cpwl = 4267 +# gamma +gamwl = cpwl / cvwl + +## FOR PATCHES 1 & 2 ## + +# density +rho0wl1 = (p01 + piwl) / ((gamwl - 1) * cvwl * T01) +rho0wl2 = (p02 + piwl) / ((gamwl - 1) * cvwl * T02) +rho0wl3 = (p03 + piwl) / ((gamwl - 1) * cvwl * T03) + +# speed of sound FOR +c_wl1 = math.sqrt(gamwl * (p01 + piwl) / rho0wl1) +c_wl2 = math.sqrt(gamwl * (p02 + piwl) / rho0wl2) +c_wl3 = math.sqrt(gamwl * (p03 + piwl) / rho0wl3) + +# part for Gases - relations from IMR +Ru = 8.3144598 # Universal gas constant (J/mol-K) + +### Vapor water ### +Rv = Ru / (18.01528e-3) # Gas constant for vapor (Ru/molecular weight) (J/kg-K) +# gamma +gamwv = 1.4 +# cp +cpwv = Rv * gamwv / (gamwv - 1) +# cv +cvwv = cpwv / gamwv +# pi infinity +piwv = 0.0e0 +# qv +qvwv = 2030000 +# qv' +qvpwv = -23400 + +## FOR PATCHES 1 & 2 ## + +# density +rho0wv1 = (p01 + piwv) / ((gamwv - 1) * cvwv * T01) +rho0wv2 = (p02 + piwv) / ((gamwv - 1) * cvwv * T02) +rho0wv3 = (p03 + piwv) / ((gamwv - 1) * cvwv * T03) + +# speed of sound +c_wv1 = math.sqrt(gamwv * (p01 + piwv) / rho0wv1) +c_wv2 = math.sqrt(gamwv * (p02 + piwv) / rho0wv2) +c_wv3 = math.sqrt(gamwv * (p03 + piwv) / rho0wv3) + +### Air ### + +Ra = Ru / (28.966e-3) # Gas constant for air (Ru/molecular weight) (J/kg-K) +gamwa = 1.4 +# cp +cpa = Ra * gamwa / (gamwa - 1) +# cv +cva = cpa / gamwa +# pi infinity +pia = 0.0e0 +# qv +qvwa = 0.0e0 +# qv' +qvpwa = 0.0e0 + +## FOR PATCHES 1 & 2 ## + +# density +rho0wa1 = (p01 + pia) / ((gamwa - 1) * cva * T01) +rho0wa2 = (p02 + pia) / ((gamwa - 1) * cva * T02) +rho0wa3 = (p03 + pia) / ((gamwa - 1) * cva * T03) + +# Speed of sound +c_a1 = math.sqrt(gamwa * (p01 + pia) / rho0wa1) +c_a2 = math.sqrt(gamwa * (p02 + pia) / rho0wa2) +c_a3 = math.sqrt(gamwa * (p03 + pia) / rho0wa3) + +### 3% polyacrylamide gel ### + +# gamma +gamwg = gamwl # 2.35 +# pi infty +pig = piwl # 1.1754E+09 +# qv +qvwg = qvwl # 0.0E0 +# qv' +qvpwg = qvpwl # 0.0E0 +# cv +cvg = cvwl +# cp +cpg = gamwg * cvg + +## FOR PATCHES 1 & 2 & 3 ## + +# density +rho0wg1 = (p01 + pig) / ((gamwg - 1) * cvg * T01) +rho0wg2 = (p02 + pig) / ((gamwg - 1) * cvg * T02) +rho0wg3 = (p03 + pig) / ((gamwg - 1) * cvg * T03) + +# Speed of sound +c_g1 = math.sqrt(gamwg * (p01 + pig) / rho0wg1) +c_g2 = math.sqrt(gamwg * (p02 + pig) / rho0wg2) +c_g3 = math.sqrt(gamwg * (p03 + pig) / rho0wg3) + +## SHOCK RELATIONS +p02Op01 = p02 / p01 + +# Mach number of the shocked region - this should agree with Min, if everything is correct +Ms = math.sqrt((gamwa + 1.0) / (2.0 * gamwa) * (p02Op01 - 1.0) * (p02 / (p02 + pia)) + 1.0) + +# shock speed +ss = Ms * c_a1 + +### volume fractions for each of the patches ### +C0 = 0.25 # vapor concentration for IMR + +# patch 1: liquid water +liq_wg = 0 +liq_wa = 1.00e-15 +liq_wv = 1.00e-15 +liq_wl = 1.00e00 - liq_wv - liq_wa - liq_wg +# water vapor +vap_wl = 1.00e-15 +vap_wv = 1 / ((1 - C0) / C0 * rho0wv2 / rho0wa2 + 1) +vap_wa = 1.00e-15 +vap_wg = 0 +vap_tot = vap_wl + vap_wv + vap_wa + vap_wg +# air +air_wl = 1.00e-15 +air_wv = vap_tot +air_wg = 0 +air_wa = 1.00e00 - air_wl - air_wv - air_wg +# bubble +bub_wl = 1e-15 +bub_wv = vap_tot +bub_wg = 0 +bub_wa = 1 - bub_wl - bub_wv - bub_wg +# gel +gel_wl = 0 +gel_wv = 0 +gel_wa = 0 +gel_wg = 1.00e00 - gel_wl - gel_wv - gel_wa + +## Elasticity +Gl = 0 +Gv = 0 +Ga = 0 +Gg = 0.57e03 + +## SIMULATION PARAMETERS + +# CFL +cfl = 0.30 + +# Bubble Initial Radius +R0 = 244.8e-06 + +# number of elements +Nx = 249 # 404 #249 +Ny = 124 # 179 #124 +Nz = 124 # 179 #124 +Nx0 = Nx + +# domain boundaries +lref = 4 * R0 +xb = -5 * R0 +xe = lref + +yb = 0.00 +ye = lref + +zb = 0.00 +ze = lref + +lenx = xe - xb +leny = ye - yb +lenz = ze - zb + +xcenl = 0.0 +ycenl = leny / 2.0 +zcenl = lenz / 2.0 + +# xdist = 6.51E-10 #2.17E-5 +# sod = xdist/R0 +sod = 1.39 +xcenb = sod * R0 # neg for bub in liq; pos bub in gel +ycenb = 0.00 +zcenb = 0.00 + +xceng = xe / 2.0 +yceng = ycenl +zceng = zcenl + +# typical cell size +dx = (xe - xb) / Nx +dy = (ye - yb) / Ny +dz = (ze - zb) / Nz +# print(dx) +PPBR_x = R0 / dx +PPBR_y = R0 / dy +PPBR_z = R0 / dz +# print(PPBR_x) +# print(PPBR_y) + +# save frequency = SF + 1 (because the initial state, 0.dat, is also saved) +SF = 100 + +# Critical time-step +tc = 0.915 * R0 * math.sqrt(rho0wl1 / p01) + +# making Nt divisible by SF +# tendA = 1.5 * tc +tend = 1.2 * tc + +# 1 - ensure NtA is sufficient to go a little beyond tendA +# NtA = int( tendA // dt + 1 ) + +# Array of saves. it is the same as Nt/Sf = t_step_save +# AS = int( NtA // SF + 1 ) + +# Nt = total number of steps. Ensure Nt > NtA (so the total tendA is covered) +# Nt = AS * SF +# Nt = int(2.5E3 * tend // tc * Nx / Nx0 + 1) +Nt = int(2.0e3 * tend // tc * Nx / Nx0 + 1) +# print(Nt) +dt = tend / Nt + +AS = int(Nt // SF) +tstart = 0 # 2184 +# Total physical time +# tend = Nt * dt + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": xb, + "x_domain%end": xe, + "y_domain%beg": yb, + "y_domain%end": ye, + "z_domain%beg": zb, + "z_domain%end": ze, + "stretch_x": "T", + "loops_x": 1, + "a_x": 4.0e0, + "x_a": -1.75 * R0 * (abs(sod) + 1), + "x_b": 5 * R0, + "stretch_y": "T", + "loops_y": 1, + "a_y": 1.0e0, + "y_a": -1.5 * R0 * abs(sod), + "y_b": 1.5 * R0 * abs(sod), + "stretch_z": "T", + "loops_z": 1, + "a_z": 1.0e0, + "z_a": -1.5 * R0 * abs(sod), + "z_b": 1.5 * R0 * abs(sod), + "cyl_coord": "F", + "m": Nx, + "n": Ny, + "p": Nz, + "dt": dt, + "t_step_start": tstart, + "t_step_stop": Nt, + "t_step_save": AS, + # Simulation Algorithm Parameters + "num_patches": 3, + "model_eqns": 3, + "num_fluids": 4, + "hypoelasticity": "F", + "hyperelasticity": "T", + "mpp_lim": "T", + "mixture_err": "T", + "relax": "T", + "relax_model": 6, + "palpha_eps": 1.0e-6, + "ptgalpha_eps": 1.0e-2, + "time_stepper": 3, + "weno_order": 3, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -6, # -2, + "bc_x%end": -6, + "bc_y%beg": -2, + "bc_y%end": -6, + "bc_z%beg": -2, + "bc_z%end": -6, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "T", + "probe_wrt": "T", + "fd_order": 1, + "num_probes": 1, + "probe(1)%x": 0.0, + "probe(1)%y": 0.0, + "probe(1)%z": 0.0, + # Patch 1: High pressured water + # Specify the cubic water background grid geometry + "patch_icpp(1)%geometry": 9, + "patch_icpp(1)%x_centroid": 20 * xcenl, + "patch_icpp(1)%y_centroid": 20 * ycenl, + "patch_icpp(1)%z_centroid": 20 * zcenl, + "patch_icpp(1)%length_x": 20 * lenx, + "patch_icpp(1)%length_y": 20 * leny, + "patch_icpp(1)%length_z": 20 * lenz, + "patch_icpp(1)%vel(1)": 0.0e00, + "patch_icpp(1)%vel(2)": 0.0e00, + "patch_icpp(1)%vel(3)": 0.0e00, + "patch_icpp(1)%pres": p01, + "patch_icpp(1)%alpha_rho(4)": liq_wl * rho0wl1, + "patch_icpp(1)%alpha_rho(2)": liq_wv * rho0wv1, + "patch_icpp(1)%alpha_rho(3)": liq_wa * rho0wa1, + "patch_icpp(1)%alpha_rho(1)": liq_wg * rho0wg1, + "patch_icpp(1)%alpha(4)": liq_wl, + "patch_icpp(1)%alpha(2)": liq_wv, + "patch_icpp(1)%alpha(3)": liq_wa, + "patch_icpp(1)%alpha(1)": liq_wg, + # Patch 2: (Vapor) Bubble + "patch_icpp(2)%geometry": 8, + "patch_icpp(2)%x_centroid": xcenb, + "patch_icpp(2)%y_centroid": ycenb, + "patch_icpp(2)%z_centroid": zcenb, + "patch_icpp(2)%radius": R0, + "patch_icpp(2)%vel(1)": 0.0e00, + "patch_icpp(2)%vel(2)": 0.0e00, + "patch_icpp(2)%vel(3)": 0.0e00, + "patch_icpp(2)%pres": p02, + "patch_icpp(2)%alpha_rho(4)": bub_wl * rho0wl2, + "patch_icpp(2)%alpha_rho(2)": bub_wv * rho0wv2, + "patch_icpp(2)%alpha_rho(3)": bub_wa * rho0wa2, + "patch_icpp(2)%alpha_rho(1)": bub_wg * rho0wg2, + "patch_icpp(2)%alpha(4)": bub_wl, + "patch_icpp(2)%alpha(2)": bub_wv, + "patch_icpp(2)%alpha(3)": bub_wa, + "patch_icpp(2)%alpha(1)": bub_wg, + "patch_icpp(2)%alter_patch(1)": "T", + # Patch 3: Gel Object + "patch_icpp(3)%geometry": 9, + "patch_icpp(3)%x_centroid": 20 * xceng, + "patch_icpp(3)%y_centroid": 20 * yceng, + "patch_icpp(3)%z_centroid": 20 * zceng, + "patch_icpp(3)%length_x": 20 * xe, + "patch_icpp(3)%length_y": 20 * leny, + "patch_icpp(3)%length_z": 20 * lenz, + "patch_icpp(3)%vel(1)": 0.0e00, + "patch_icpp(3)%vel(2)": 0.0e00, + "patch_icpp(3)%vel(3)": 0.0e00, + "patch_icpp(3)%pres": p03, + "patch_icpp(3)%alpha_rho(4)": gel_wl * rho0wl3, + "patch_icpp(3)%alpha_rho(2)": gel_wv * rho0wv3, + "patch_icpp(3)%alpha_rho(3)": gel_wa * rho0wa3, + "patch_icpp(3)%alpha_rho(1)": gel_wg * rho0wg3, + "patch_icpp(3)%alpha(4)": gel_wl, + "patch_icpp(3)%alpha(2)": gel_wv, + "patch_icpp(3)%alpha(3)": gel_wa, + "patch_icpp(3)%alpha(1)": gel_wg, + "patch_icpp(3)%alter_patch(1)": "T", + # Fluids Physical Parameters + "fluid_pp(4)%gamma": 1.0e00 / (gamwl - 1), + "fluid_pp(4)%pi_inf": gamwl * piwl / (gamwl - 1), + "fluid_pp(4)%cv": cvwl, + "fluid_pp(4)%qv": qvwl, + "fluid_pp(4)%qvp": qvpwl, + "fluid_pp(4)%G": Gl, + "fluid_pp(2)%gamma": 1.0e00 / (gamwv - 1), + "fluid_pp(2)%pi_inf": gamwv * piwv / (gamwv - 1), + "fluid_pp(2)%cv": cvwv, + "fluid_pp(2)%qv": qvwv, + "fluid_pp(2)%qvp": qvpwv, + "fluid_pp(2)%G": Gv, + "fluid_pp(3)%gamma": 1.0e00 / (gamwa - 1), + "fluid_pp(3)%pi_inf": gamwa * pia / (gamwa - 1), + "fluid_pp(3)%cv": cva, + "fluid_pp(3)%qv": qvwa, + "fluid_pp(3)%qvp": qvpwa, + "fluid_pp(3)%G": Ga, + "fluid_pp(1)%gamma": 1.0e00 / (gamwg - 1), + "fluid_pp(1)%pi_inf": gamwg * pig / (gamwg - 1), + "fluid_pp(1)%cv": cvg, + "fluid_pp(1)%qv": qvwg, + "fluid_pp(1)%qvp": qvpwg, + "fluid_pp(1)%G": Gg, + } + ) +) diff --git a/examples/3D_hyper_bubinwater/case.py b/examples/3D_hyper_bubinwater/case.py new file mode 100755 index 0000000000..4e01000605 --- /dev/null +++ b/examples/3D_hyper_bubinwater/case.py @@ -0,0 +1,414 @@ +#!/usr/bin/env python3 +import math, json + +## 1 FOR BACKGROUND, 2 FOR BUBBLE, 3 FOR GEL +# Pressure [Pa] +p01 = 5e6 +p02 = 3550 +p03 = p01 + +# Temperature [K] +T01 = 298.15 +T02 = 298.15 +T03 = T01 + +#### FLUID PROPERTIES #### + +### liquid water ### +# pi infty +piwl = 1.0e09 +# qv +qvwl = -1167000 +# qv' +qvpwl = 0.0e0 +# cv +cvwl = 1816 +# cp +cpwl = 4267 +# gamma +gamwl = cpwl / cvwl + +## FOR PATCHES 1 & 2 ## + +# density +rho0wl1 = (p01 + piwl) / ((gamwl - 1) * cvwl * T01) +rho0wl2 = (p02 + piwl) / ((gamwl - 1) * cvwl * T02) +rho0wl3 = (p03 + piwl) / ((gamwl - 1) * cvwl * T03) + +# speed of sound FOR +c_wl1 = math.sqrt(gamwl * (p01 + piwl) / rho0wl1) +c_wl2 = math.sqrt(gamwl * (p02 + piwl) / rho0wl2) +c_wl3 = math.sqrt(gamwl * (p03 + piwl) / rho0wl3) + +# part for Gases - relations from IMR +Ru = 8.3144598 # Universal gas constant (J/mol-K) + +### Vapor water ### +Rv = Ru / (18.01528e-3) # Gas constant for vapor (Ru/molecular weight) (J/kg-K) +# gamma +gamwv = 1.4 +# cp +cpwv = Rv * gamwv / (gamwv - 1) +# cv +cvwv = cpwv / gamwv +# pi infinity +piwv = 0.0e0 +# qv +qvwv = 2030000 +# qv' +qvpwv = -23400 + +## FOR PATCHES 1 & 2 ## + +# density +rho0wv1 = (p01 + piwv) / ((gamwv - 1) * cvwv * T01) +rho0wv2 = (p02 + piwv) / ((gamwv - 1) * cvwv * T02) +rho0wv3 = (p03 + piwv) / ((gamwv - 1) * cvwv * T03) + +# speed of sound +c_wv1 = math.sqrt(gamwv * (p01 + piwv) / rho0wv1) +c_wv2 = math.sqrt(gamwv * (p02 + piwv) / rho0wv2) +c_wv3 = math.sqrt(gamwv * (p03 + piwv) / rho0wv3) + +### Air ### + +Ra = Ru / (28.966e-3) # Gas constant for air (Ru/molecular weight) (J/kg-K) +gamwa = 1.4 +# cp +cpa = Ra * gamwa / (gamwa - 1) +# cv +cva = cpa / gamwa +# pi infinity +pia = 0.0e0 +# qv +qvwa = 0.0e0 +# qv' +qvpwa = 0.0e0 + +## FOR PATCHES 1 & 2 ## + +# density +rho0wa1 = (p01 + pia) / ((gamwa - 1) * cva * T01) +rho0wa2 = (p02 + pia) / ((gamwa - 1) * cva * T02) +rho0wa3 = (p03 + pia) / ((gamwa - 1) * cva * T03) + +# Speed of sound +c_a1 = math.sqrt(gamwa * (p01 + pia) / rho0wa1) +c_a2 = math.sqrt(gamwa * (p02 + pia) / rho0wa2) +c_a3 = math.sqrt(gamwa * (p03 + pia) / rho0wa3) + +### 3% polyacrylamide gel ### + +# gamma +gamwg = gamwl # 2.35 +# pi infty +pig = piwl # 1.1754E+09 +# qv +qvwg = qvwl # 0.0E0 +# qv' +qvpwg = qvpwl # 0.0E0 +# cv +cvg = cvwl +# cp +cpg = gamwg * cvg + +## FOR PATCHES 1 & 2 & 3 ## + +# density +rho0wg1 = (p01 + pig) / ((gamwg - 1) * cvg * T01) +rho0wg2 = (p02 + pig) / ((gamwg - 1) * cvg * T02) +rho0wg3 = (p03 + pig) / ((gamwg - 1) * cvg * T03) + +# Speed of sound +c_g1 = math.sqrt(gamwg * (p01 + pig) / rho0wg1) +c_g2 = math.sqrt(gamwg * (p02 + pig) / rho0wg2) +c_g3 = math.sqrt(gamwg * (p03 + pig) / rho0wg3) + +## SHOCK RELATIONS +p02Op01 = p02 / p01 + +# Mach number of the shocked region - this should agree with Min, if everything is correct +Ms = math.sqrt((gamwa + 1.0) / (2.0 * gamwa) * (p02Op01 - 1.0) * (p02 / (p02 + pia)) + 1.0) + +# shock speed +ss = Ms * c_a1 + +### volume fractions for each of the patches ### +C0 = 0.25 # vapor concentration for IMR + +# patch 1: liquid water +liq_wg = 0 +liq_wa = 1.00e-15 +liq_wv = 1.00e-15 +liq_wl = 1.00e00 - liq_wv - liq_wa - liq_wg +# water vapor +vap_wl = 1.00e-15 +vap_wv = 1 / ((1 - C0) / C0 * rho0wv2 / rho0wa2 + 1) +vap_wa = 1.00e-15 +vap_wg = 0 +vap_tot = vap_wl + vap_wv + vap_wa + vap_wg +# air +air_wl = 1.00e-15 +air_wv = vap_tot +air_wg = 0 +air_wa = 1.00e00 - air_wl - air_wv - air_wg +# bubble +bub_wl = 1e-15 +bub_wv = vap_tot +bub_wg = 0 +bub_wa = 1 - bub_wl - bub_wv - bub_wg +# gel +gel_wl = 0 +gel_wv = 0 +gel_wa = 0 +gel_wg = 1.00e00 - gel_wl - gel_wv - gel_wa + +## Elasticity +Gl = 0 +Gv = 0 +Ga = 0 +Gg = 0.57e03 + +## SIMULATION PARAMETERS + +# CFL +cfl = 0.50 + +# Bubble Initial Radius +R0 = 230.4e-06 + +# number of elements +Nx = 249 # 404 #249 +Ny = 124 # 179 #124 +Nz = 124 # 179 #124 +Nx0 = Nx + +# domain boundaries +lref = 4 * R0 +xb = -5 * R0 +xe = lref + +yb = 0.00 +ye = lref + +zb = 0.00 +ze = lref + +lenx = xe - xb +leny = ye - yb +lenz = ze - zb + +xcenl = 0.0 +ycenl = leny / 2.0 +zcenl = lenz / 2.0 + +# xdist = 6.51E-10 #2.17E-5 +# sod = xdist/R0 +sod = 2.17 +xcenb = -sod * R0 # -sod #-sod*R0 +ycenb = 0.00 +zcenb = 0.00 + +xceng = xe / 2.0 +yceng = ycenl +zceng = zcenl + +# typical cell size +dx = (xe - xb) / Nx +dy = (ye - yb) / Ny +dz = (ze - zb) / Nz +# print(dx) +# time step + +# save frequency = SF + 1 (because the initial state, 0.dat, is also saved) +SF = 100 + +# Critical time-step +tc = 0.915 * R0 * math.sqrt(rho0wl1 / p01) + +# making Nt divisible by SF +# tendA = 1.5 * tc +tend = 1.2 * tc + +# 1 - ensure NtA is sufficient to go a little beyond tendA +# NtA = int( tendA // dt + 1 ) + +# Array of saves. it is the same as Nt/Sf = t_step_save +# AS = int( NtA // SF + 1 ) + +# Nt = total number of steps. Ensure Nt > NtA (so the total tendA is covered) +# Nt = AS * SF +Nt = int(2e3 * tend // tc * Nx / Nx0 + 1) +# print(Nt) +dt = tend / Nt + +AS = int(Nt // SF) +tstart = 0 # 2304 +# Total physical time +# tend = Nt * dt + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": xb, + "x_domain%end": xe, + "y_domain%beg": yb, + "y_domain%end": ye, + "z_domain%beg": zb, + "z_domain%end": ze, + "stretch_x": "T", + "loops_x": 1, + "a_x": 4.0e0, + "x_a": -1.75 * R0 * (abs(sod) + 1), + "x_b": 5 * R0, + "stretch_y": "T", + "loops_y": 1, + "a_y": 1.0e0, + "y_a": -1.5 * R0 * abs(sod), + "y_b": 1.5 * R0 * abs(sod), + "stretch_z": "T", + "loops_z": 1, + "a_z": 1.0e0, + "z_a": -1.5 * R0 * abs(sod), + "z_b": 1.5 * R0 * abs(sod), + "cyl_coord": "F", + "m": Nx, + "n": Ny, + "p": Nz, + "dt": dt, + "t_step_start": tstart, + "t_step_stop": Nt, + "t_step_save": AS, + # Simulation Algorithm Parameters + "num_patches": 3, + "model_eqns": 3, + "num_fluids": 4, + "hypoelasticity": "F", + "hyperelasticity": "T", + "mpp_lim": "T", + "mixture_err": "T", + "relax": "T", + "relax_model": 6, + "palpha_eps": 1.0e-6, + "ptgalpha_eps": 1.0e-2, + "time_stepper": 3, + "weno_order": 3, + "weno_eps": 1.0e-16, + "weno_Re_flux": "F", + "weno_avg": "F", + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + "riemann_solver": 2, + "wave_speeds": 1, + "avg_state": 2, + "bc_x%beg": -6, # -2, + "bc_x%end": -6, + "bc_y%beg": -2, + "bc_y%end": -6, + "bc_z%beg": -2, + "bc_z%end": -6, + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "T", + "probe_wrt": "T", + "fd_order": 1, + "num_probes": 1, + "probe(1)%x": 0.0, + "probe(1)%y": 0.0, + "probe(1)%z": 0.0, + # Patch 1: High pressured water + # Specify the cubic water background grid geometry + "patch_icpp(1)%geometry": 9, + "patch_icpp(1)%x_centroid": 20 * xcenl, + "patch_icpp(1)%y_centroid": 20 * ycenl, + "patch_icpp(1)%z_centroid": 20 * zcenl, + "patch_icpp(1)%length_x": 20 * lenx, + "patch_icpp(1)%length_y": 20 * leny, + "patch_icpp(1)%length_z": 20 * lenz, + "patch_icpp(1)%vel(1)": 0.0e00, + "patch_icpp(1)%vel(2)": 0.0e00, + "patch_icpp(1)%vel(3)": 0.0e00, + "patch_icpp(1)%pres": p01, + "patch_icpp(1)%alpha_rho(1)": liq_wl * rho0wl1, + "patch_icpp(1)%alpha_rho(2)": liq_wv * rho0wv1, + "patch_icpp(1)%alpha_rho(3)": liq_wa * rho0wa1, + "patch_icpp(1)%alpha_rho(4)": liq_wg * rho0wg1, + "patch_icpp(1)%alpha(1)": liq_wl, + "patch_icpp(1)%alpha(2)": liq_wv, + "patch_icpp(1)%alpha(3)": liq_wa, + "patch_icpp(1)%alpha(4)": liq_wg, + # Patch 2: (Vapor) Bubble + "patch_icpp(2)%geometry": 8, + "patch_icpp(2)%x_centroid": xcenb, + "patch_icpp(2)%y_centroid": ycenb, + "patch_icpp(2)%z_centroid": zcenb, + "patch_icpp(2)%radius": R0, + "patch_icpp(2)%vel(1)": 0.0e00, + "patch_icpp(2)%vel(2)": 0.0e00, + "patch_icpp(2)%vel(3)": 0.0e00, + "patch_icpp(2)%pres": p02, + "patch_icpp(2)%alpha_rho(1)": bub_wl * rho0wl2, + "patch_icpp(2)%alpha_rho(2)": bub_wv * rho0wv2, + "patch_icpp(2)%alpha_rho(3)": bub_wa * rho0wa2, + "patch_icpp(2)%alpha_rho(4)": bub_wg * rho0wg2, + "patch_icpp(2)%alpha(1)": bub_wl, + "patch_icpp(2)%alpha(2)": bub_wv, + "patch_icpp(2)%alpha(3)": bub_wa, + "patch_icpp(2)%alpha(4)": bub_wg, + "patch_icpp(2)%alter_patch(1)": "T", + # Patch 3: Gel Object + "patch_icpp(3)%geometry": 9, + "patch_icpp(3)%x_centroid": 20 * xceng, + "patch_icpp(3)%y_centroid": 20 * yceng, + "patch_icpp(3)%z_centroid": 20 * zceng, + "patch_icpp(3)%length_x": 20 * xe, + "patch_icpp(3)%length_y": 20 * leny, + "patch_icpp(3)%length_z": 20 * lenz, + "patch_icpp(3)%vel(1)": 0.0e00, + "patch_icpp(3)%vel(2)": 0.0e00, + "patch_icpp(3)%vel(3)": 0.0e00, + "patch_icpp(3)%pres": p03, + "patch_icpp(3)%alpha_rho(1)": gel_wl * rho0wl3, + "patch_icpp(3)%alpha_rho(2)": gel_wv * rho0wv3, + "patch_icpp(3)%alpha_rho(3)": gel_wa * rho0wa3, + "patch_icpp(3)%alpha_rho(4)": gel_wg * rho0wg3, + "patch_icpp(3)%alpha(1)": gel_wl, + "patch_icpp(3)%alpha(2)": gel_wv, + "patch_icpp(3)%alpha(3)": gel_wa, + "patch_icpp(3)%alpha(4)": gel_wg, + "patch_icpp(3)%alter_patch(1)": "T", + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (gamwl - 1), + "fluid_pp(1)%pi_inf": gamwl * piwl / (gamwl - 1), + "fluid_pp(1)%cv": cvwl, + "fluid_pp(1)%qv": qvwl, + "fluid_pp(1)%qvp": qvpwl, + "fluid_pp(1)%G": Gl, + "fluid_pp(2)%gamma": 1.0e00 / (gamwv - 1), + "fluid_pp(2)%pi_inf": gamwv * piwv / (gamwv - 1), + "fluid_pp(2)%cv": cvwv, + "fluid_pp(2)%qv": qvwv, + "fluid_pp(2)%qvp": qvpwv, + "fluid_pp(2)%G": Gv, + "fluid_pp(3)%gamma": 1.0e00 / (gamwa - 1), + "fluid_pp(3)%pi_inf": gamwa * pia / (gamwa - 1), + "fluid_pp(3)%cv": cva, + "fluid_pp(3)%qv": qvwa, + "fluid_pp(3)%qvp": qvpwa, + "fluid_pp(3)%G": Ga, + "fluid_pp(4)%gamma": 1.0e00 / (gamwg - 1), + "fluid_pp(4)%pi_inf": gamwg * pig / (gamwg - 1), + "fluid_pp(4)%cv": cvg, + "fluid_pp(4)%qv": qvwg, + "fluid_pp(4)%qvp": qvpwg, + "fluid_pp(4)%G": Gg, + } + ) +) diff --git a/src/common/m_boundary_conditions.fpp b/src/common/m_boundary_conditions.fpp index 762ff6a44f..87be1d10ad 100644 --- a/src/common/m_boundary_conditions.fpp +++ b/src/common/m_boundary_conditions.fpp @@ -27,6 +27,7 @@ contains !> The purpose of this procedure is to populate the buffers !! of the primitive variables, depending on the selected !! boundary conditions. + !! @param q_prim_vf Primitive variable subroutine s_populate_variables_buffers(q_prim_vf, pb, mv) type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf @@ -218,6 +219,7 @@ contains real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(inout) :: pb, mv integer, intent(in) :: bc_dir, bc_loc integer :: j, k, l, q, i + real(wp) :: bc_sum !< x-direction if (bc_dir == 1) then !< x-direction @@ -235,7 +237,27 @@ contains end do end do end do - +#ifdef MFC_SIMULATION + if (hyperelasticity) then + !$acc parallel loop collapse(4) gang vector default(present), private(bc_sum) + do j = 1, buff_size + do l = 0, p + do k = 0, n + do i = xibeg, xiend + bc_sum = 0._wp + !$acc loop seq + do q = 1, j + bc_sum = bc_sum - dx(-q) + end do + q_prim_vf(i)%sf(-j, k, l) = & + q_prim_vf(i)%sf(0, k, l) - bc_sum + end do + end do + end do + end do + !$acc end parallel loop + end if +#endif else !< bc_x%end !$acc parallel loop collapse(4) gang vector default(present) @@ -249,7 +271,27 @@ contains end do end do end do - +#ifdef MFC_SIMULATION + if (hyperelasticity) then + !$acc parallel loop collapse(4) gang vector default(present), private(bc_sum) + do j = 1, buff_size + do l = 0, p + do k = 0, n + do i = xibeg, xiend + bc_sum = 0._wp + !$acc loop seq + do q = 1, j + bc_sum = bc_sum + dx(m + q) + end do + q_prim_vf(i)%sf(m + j, k, l) = & + q_prim_vf(i)%sf(m, k, l) + bc_sum + end do + end do + end do + end do + !$acc end parallel loop + end if +#endif end if !< y-direction @@ -268,7 +310,27 @@ contains end do end do end do - +#ifdef MFC_SIMULATION + if (hyperelasticity) then + !$acc parallel loop collapse(4) gang vector default(present), private(bc_sum) + do j = 1, buff_size + do l = -buff_size, m + buff_size + do k = 0, p + do i = xibeg, xiend + bc_sum = 0._wp + !$acc loop seq + do q = 1, j + bc_sum = bc_sum - dy(-q) + end do + q_prim_vf(i)%sf(l, -j, k) = & + q_prim_vf(i)%sf(l, 0, k) - bc_sum + end do + end do + end do + end do + !$acc end parallel loop + end if +#endif else !< bc_y%end !$acc parallel loop collapse(4) gang vector default(present) @@ -282,7 +344,27 @@ contains end do end do end do - +#ifdef MFC_SIMULATION + if (hyperelasticity) then + !$acc parallel loop collapse(4) gang vector default(present), private(bc_sum) + do j = 1, buff_size + do l = -buff_size, m + buff_size + do k = 0, p + do i = xibeg, xiend + bc_sum = 0._wp + !$acc loop seq + do q = 1, j + bc_sum = bc_sum + dy(n + q) + end do + q_prim_vf(i)%sf(l, n + j, k) = & + q_prim_vf(i)%sf(l, n, k) + bc_sum + end do + end do + end do + end do + !$acc end parallel loop + end if +#endif end if !< z-direction @@ -301,9 +383,28 @@ contains end do end do end do - +#ifdef MFC_SIMULATION + if (hyperelasticity) then + !$acc parallel loop collapse(4) gang vector default(present), private(bc_sum) + do j = 1, buff_size + do l = -buff_size, n + buff_size + do k = -buff_size, m + buff_size + do i = xibeg, xiend + bc_sum = 0._wp + !$acc loop seq + do q = 1, j + bc_sum = bc_sum - dx(-q) + end do + q_prim_vf(i)%sf(k, l, -j) = & + q_prim_vf(i)%sf(k, l, 0) - bc_sum + end do + end do + end do + end do + !$acc end parallel loop + end if +#endif else !< bc_z%end - !$acc parallel loop collapse(4) gang vector default(present) do i = 1, sys_size do j = 1, buff_size @@ -315,7 +416,27 @@ contains end do end do end do - +#ifdef MFC_SIMULATION + if (hyperelasticity) then + !$acc parallel loop collapse(4) gang vector default(present), private(bc_sum) + do j = 1, buff_size + do l = -buff_size, n + buff_size + do k = -buff_size, m + buff_size + do i = xibeg, xiend + bc_sum = 0._wp + !$acc loop seq + do q = 1, j + bc_sum = bc_sum + dz(p + q) + end do + q_prim_vf(i)%sf(k, l, p + j) = & + q_prim_vf(i)%sf(k, l, p) + bc_sum + end do + end do + end do + end do + !$acc end parallel loop + end if +#endif end if end if @@ -361,11 +482,6 @@ contains end do end if - if (hyperelasticity) then - q_prim_vf(xibeg)%sf(-j, k, l) = & - -q_prim_vf(xibeg)%sf(j - 1, k, l) - end if - end do end do end do @@ -418,11 +534,6 @@ contains end do end if - if (hyperelasticity) then - q_prim_vf(xibeg)%sf(m + j, k, l) = & - -q_prim_vf(xibeg)%sf(m - (j - 1), k, l) - end if - end do end do end do @@ -478,10 +589,6 @@ contains end do end if - if (hyperelasticity) then - q_prim_vf(xibeg + 1)%sf(l, -j, k) = & - -q_prim_vf(xibeg + 1)%sf(l, j - 1, k) - end if end do end do end do @@ -532,10 +639,6 @@ contains end do end if - if (hyperelasticity) then - q_prim_vf(xibeg + 1)%sf(l, n + j, k) = & - -q_prim_vf(xibeg + 1)%sf(l, n - (j - 1), k) - end if end do end do end do @@ -591,10 +694,6 @@ contains end do end if - if (hyperelasticity) then - q_prim_vf(xiend)%sf(k, l, -j) = & - -q_prim_vf(xiend)%sf(k, l, j - 1) - end if end do end do end do @@ -645,10 +744,6 @@ contains end do end if - if (hyperelasticity) then - q_prim_vf(xiend)%sf(k, l, p + j) = & - -q_prim_vf(xiend)%sf(k, l, p - (j - 1)) - end if end do end do end do diff --git a/src/common/m_checker_common.fpp b/src/common/m_checker_common.fpp index 483de492dd..2281531852 100644 --- a/src/common/m_checker_common.fpp +++ b/src/common/m_checker_common.fpp @@ -39,7 +39,6 @@ contains call s_check_inputs_bubbles_euler call s_check_inputs_qbmm_and_polydisperse call s_check_inputs_adv_n - call s_check_inputs_hypoelasticity call s_check_inputs_phase_change call s_check_inputs_ibm #endif @@ -50,6 +49,7 @@ contains call s_check_inputs_weno call s_check_inputs_bc call s_check_inputs_stiffened_eos + call s_check_inputs_elasticity call s_check_inputs_surface_tension call s_check_inputs_moving_bc call s_check_inputs_mhd @@ -133,25 +133,6 @@ contains @:PROHIBIT(adv_n .and. qbmm) end subroutine s_check_inputs_adv_n - !> Checks constraints on the hypoelasticity parameters. - !! Called by s_check_inputs_common for pre-processing and simulation - subroutine s_check_inputs_hypoelasticity - @:PROHIBIT(hypoelasticity .and. model_eqns /= 2) -#ifdef MFC_SIMULATION - @:PROHIBIT(elasticity .and. fd_order /= 4) -#endif - end subroutine s_check_inputs_hypoelasticity - - !> Checks constraints on the hyperelasticity parameters. - !! Called by s_check_inputs_common for pre-processing and simulation - subroutine s_check_inputs_hyperelasticity - @:PROHIBIT(hyperelasticity .and. model_eqns == 1) - @:PROHIBIT(hyperelasticity .and. model_eqns > 3) -#ifdef MFC_SIMULATION - @:PROHIBIT(elasticity .and. fd_order /= 4) -#endif - end subroutine s_check_inputs_hyperelasticity - !> Checks constraints on the phase change parameters. !! Called by s_check_inputs_common for pre-processing and simulation subroutine s_check_inputs_phase_change @@ -178,6 +159,17 @@ contains #endif + !> Checks constraints on the elasticity parameters. + !! Called by s_check_inputs_common for all three stages + subroutine s_check_inputs_elasticity + @:PROHIBIT(elasticity .and. .not. (model_eqns == 2 .or. model_eqns == 3)) + #:for X in ['x', 'y', 'z'] + #:for BOUND in ['beg', 'end'] + @:PROHIBIT(hyperelasticity .and. bc_${X}$%${BOUND}$ /= dflt_int .and. (bc_${X}$%${BOUND}$ < -3)) + #:endfor + #:endfor + end subroutine s_check_inputs_elasticity + !> Checks constraints on dimensionality and the number of cells for the grid. !! Called by s_check_inputs_common for all three stages subroutine s_check_inputs_simulation_domain diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 3e6b70d6b7..d0bd2cd0f0 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -1028,7 +1028,7 @@ contains if (model_eqns /= 4) then qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) & /rho_K - dyn_pres_K = dyn_pres_K + 5e-1_wp*qK_cons_vf(i)%sf(j, k, l) & + dyn_pres_K = dyn_pres_K + 0.5_wp*qK_cons_vf(i)%sf(j, k, l) & *qK_prim_vf(i)%sf(j, k, l) else qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) & @@ -1140,6 +1140,16 @@ contains end do end if +#ifdef MFC_POST_PROCESS + if (hyperelasticity .and. p /= 0) then + ! to save von Mises stress instead of elastic internal energy + qK_prim_vf(xiend + 1)%sf(j, k, l) = sqrt((3_wp/2_wp)*(qK_prim_vf(strxb)%sf(j, k, l)**2_wp + & + 2_wp*qK_prim_vf(strxb + 1)%sf(j, k, l)**2_wp + qK_prim_vf(strxb + 2)%sf(j, k, l)**2_wp + & + 2_wp*qK_prim_vf(strxb + 3)%sf(j, k, l)**2_wp + 2_wp*qK_prim_vf(strxb + 4)%sf(j, k, l)**2_wp + & + qK_prim_vf(strxe)%sf(j, k, l)**2_wp)) + end if +#endif + !$acc loop seq do i = advxb, advxe qK_prim_vf(i)%sf(j, k, l) = qK_cons_vf(i)%sf(j, k, l) @@ -1512,7 +1522,7 @@ contains ! Computing the energy from the pressure E_K = gamma_K*pres_K + pi_inf_K & - + 5e-1_wp*rho_K*vel_K_sum + qv_K + + 0.5_wp*rho_K*vel_K_sum + qv_K ! mass flux, this should be \alpha_i \rho_i u_i !$acc loop seq @@ -1640,7 +1650,7 @@ contains (rho*(1._wp - adv(num_fluids))) end if else - c = ((H - 5e-1*vel_sum)/gamma) + c = ((H - 0.5_wp*vel_sum)/gamma) end if if (mixture_err .and. c < 0._wp) then diff --git a/src/post_process/m_checker.fpp b/src/post_process/m_checker.fpp index cfe8d8ad27..8d2f23b664 100644 --- a/src/post_process/m_checker.fpp +++ b/src/post_process/m_checker.fpp @@ -33,6 +33,7 @@ contains call s_check_inputs_volume_fraction call s_check_inputs_vorticity call s_check_inputs_schlieren + call s_check_inputs_elasticity call s_check_inputs_surface_tension call s_check_inputs_no_flow_variables @@ -129,6 +130,14 @@ contains end do end subroutine s_check_inputs_schlieren + !> Checks constraints on elasticity parameters + subroutine s_check_inputs_elasticity + @:PROHIBIT(.not. (hypoelasticity .or. hyperelasticity) .and. tau_wrt) + ! Note: 'elasticity' variable isn't initialized yet; use (hypoelasticity .or. hyperelasticity) instead + @:PROHIBIT(.not. hyperelasticity .and. kymograph) + @:PROHIBIT(format == 2 .and. kymograph, 'Binary output format does not support kymograph') + end subroutine + !> Checks constraints on surface tension parameters (cf_wrt and sigma) subroutine s_check_inputs_surface_tension @:PROHIBIT(cf_wrt .and. .not. surface_tension, & @@ -138,7 +147,7 @@ contains !> Checks constraints on the absence of flow variables subroutine s_check_inputs_no_flow_variables @:PROHIBIT(.not. any([ & - (/rho_wrt, E_wrt, pres_wrt, gamma_wrt, heat_ratio_wrt, pi_inf_wrt, & + (/rho_wrt, E_wrt, pres_wrt, tau_wrt, gamma_wrt, heat_ratio_wrt, pi_inf_wrt, & pres_inf_wrt, cons_vars_wrt, prim_vars_wrt, c_wrt, schlieren_wrt/), & alpha_rho_wrt, mom_wrt, vel_wrt, flux_wrt, alpha_wrt, omega_wrt]), & "None of the flow variables have been selected for post-process. Exiting.") diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp index 41008ed7e6..e658fc14e5 100644 --- a/src/post_process/m_data_output.fpp +++ b/src/post_process/m_data_output.fpp @@ -28,19 +28,21 @@ module m_data_output s_open_formatted_database_file, & s_open_intf_data_file, & s_open_energy_data_file, & + s_open_kymo_data_file, & s_write_grid_to_formatted_database_file, & s_write_variable_to_formatted_database_file, & s_write_lag_bubbles_results, & s_write_intf_data_file, & s_write_energy_data_file, & + s_write_kymo_data_file, & s_close_formatted_database_file, & s_close_intf_data_file, & s_close_energy_data_file, & + s_close_kymo_data_file, & s_finalize_data_output_module ! Including the Silo Fortran interface library that features the subroutines ! and parameters that are required to write in the Silo-HDF5 database format - ! INCLUDE 'silo.inc' include 'silo_f9x.inc' ! Generic storage for flow variable(s) that are to be written to formatted @@ -344,7 +346,7 @@ contains if (pres_wrt .or. prim_vars_wrt) dbvars = dbvars + 1 ! Elastic stresses - if (hypoelasticity) dbvars = dbvars + (num_dims*(num_dims + 1))/2 + if (tau_wrt .or. prim_vars_wrt) dbvars = dbvars + (num_dims*(num_dims + 1))/2 ! Magnetic field if (mhd) then @@ -614,6 +616,25 @@ contains end subroutine s_open_energy_data_file + subroutine s_open_kymo_data_file() + ! Time-step that is currently being post-processed + + ! Relative path to a file in the case directory + character(LEN=path_len + 3*name_len) :: file_path + + ! Kymo information is in binary database format + ! Generates relative path to database, opened for current time-step + write (file_path, '(A)') '/kymo_data.dat' + file_path = trim(case_dir)//trim(file_path) + + ! Opening the simulation data file + open (251, FILE=trim(file_path), & + FORM='formatted', & + POSITION='append', & + STATUS='unknown') + + end subroutine s_open_kymo_data_file + subroutine s_write_grid_to_formatted_database_file(t_step) ! Description: The general objective of this subroutine is to write the ! necessary grid data to the formatted database file, for @@ -1381,6 +1402,42 @@ contains end subroutine s_write_energy_data_file + subroutine s_write_kymo_data_file(q_prim_vf) + type(scalar_field), dimension(sys_size), intent(IN) :: q_prim_vf + integer :: j, k, l, t !< Generic loop iterators + real(wp) :: vonMises_d, vonMises_h1 !< selected planes for kymograph comparison + real(wp) :: vonMises_h2, vonMises_h3 ! null() -! if (sim_data .and. proc_rank == 0) then -! call s_close_intf_data_file() -! call s_close_energy_data_file() -! end if - ! Deallocation procedures for the modules call s_finalize_data_output_module() call s_finalize_derived_variables_module() diff --git a/src/pre_process/include/3dHardcodedIC.fpp b/src/pre_process/include/3dHardcodedIC.fpp index e0018598a7..a59f7b3c68 100644 --- a/src/pre_process/include/3dHardcodedIC.fpp +++ b/src/pre_process/include/3dHardcodedIC.fpp @@ -4,6 +4,10 @@ real(wp) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph real(wp) :: eps + real(wp) :: rcoord, theta, phi, xi_sph, x_bcen, y_bcen, z_bcen, Rinit + real(wp) :: x_ccs, y_ccs, z_ccs + real(wp), dimension(num_dims) :: xi_cart + integer :: l eps = 1e-9_wp #:enddef @@ -44,8 +48,8 @@ end if case (301) ! (3D lung geometry in X direction, |sin(*)+sin(*)|) - h = 0.0_wp - lam = 1.0_wp + h = 0._wp + lam = 1._wp amp = patch_icpp(patch_id)%a(2) intH = amp*abs((sin(2*pi*y_cc(j)/lam - pi/2) + sin(2*pi*z_cc(k)/lam - pi/2)) + h) if (x_cc(i) > intH) then @@ -56,6 +60,60 @@ q_prim_vf(advxe)%sf(i, j, k) = patch_icpp(1)%alpha(2) end if + case (302) ! (3D lung geometry in X direction - axisym, with smoothing) + lam = 200.e-06_wp + amp = patch_icpp(patch_id)%a(2) + h = 0.125_wp*amp + + intH = amp/2._wp*(sin(2._wp*pi*y_cc(j)/lam + pi/2._wp) + sin(2._wp*pi*z_cc(k)/lam + pi/2._wp)) + + alph = patch_icpp(2)%alpha(1) + (patch_icpp(1)%alpha(1) - patch_icpp(2)%alpha(1))/(h)*(x_cc(i) - (intH - h/2._wp)) + + if (x_cc(i) > intH + h/2) then + + q_prim_vf(advxb)%sf(i, j, k) = patch_icpp(1)%alpha(1) + q_prim_vf(advxe)%sf(i, j, k) = patch_icpp(1)%alpha(2) + q_prim_vf(contxb)%sf(i, j, k) = patch_icpp(1)%alpha_rho(1) + q_prim_vf(contxe)%sf(i, j, k) = patch_icpp(1)%alpha_rho(2) + q_prim_vf(E_idx)%sf(i, j, k) = patch_icpp(1)%pres + + else if ((x_cc(i) <= intH + h/2) .and. (x_cc(i) >= intH - h/2._wp)) then + + q_prim_vf(advxb)%sf(i, j, k) = alph !0.5 + q_prim_vf(advxe)%sf(i, j, k) = 1._wp - alph !0.5 + q_prim_vf(contxb)%sf(i, j, k) = patch_icpp(1)%alpha_rho(1)/patch_icpp(1)%alpha(1)*alph!0.5 + q_prim_vf(contxe)%sf(i, j, k) = patch_icpp(2)%alpha_rho(2)/patch_icpp(2)%alpha(2)*(1 - alph)!0.5 + q_prim_vf(E_idx)%sf(i, j, k) = patch_icpp(1)%pres + + end if + + case (303) ! pre_stress for hyperelasticity, bubble in material + + R0ref = 30e-6_wp ! equilibrium radius + Rinit = patch_icpp(3)%radius ! initial radius + x_bcen = patch_icpp(3)%x_centroid + y_bcen = patch_icpp(3)%y_centroid + z_bcen = patch_icpp(3)%z_centroid + x_ccs = x_cc(i) - x_bcen + y_ccs = y_cc(j) - y_bcen + z_ccs = z_cc(k) - z_bcen + rcoord = sqrt(x_ccs**2._wp + y_ccs**2._wp + z_ccs**2._wp) + phi = atan2(y_ccs, x_ccs) + theta = atan2(sqrt(x_ccs**2._wp + y_ccs**2._wp), z_ccs) + !spherical coord, assuming Rmax=1 + xi_sph = (rcoord**3._wp - R0ref**3._wp + Rinit**3._wp)**(1._wp/3._wp) + xi_cart(1) = xi_sph*sin(theta)*cos(phi) + xi_cart(2) = xi_sph*sin(theta)*sin(phi) + xi_cart(3) = xi_sph*cos(theta) + ! shift back + xi_cart(1) = xi_cart(1) + x_bcen + xi_cart(2) = xi_cart(2) + y_bcen + xi_cart(3) = xi_cart(3) + z_bcen + ! assigning the reference map to the q_prim vector field + do l = 1, 3 + q_prim_vf(l + xibeg - 1)%sf(i, j, k) = xi_cart(l) + end do + ! Put your variable assignments here case default call s_int_to_str(patch_id, iStr) diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index ba5ab6608d..b45145dcbd 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -508,21 +508,13 @@ contains end if ! Elastic Shear Stress - if (hyperelasticity) then - - if (pre_stress) then ! pre stressed initial condition in spatial domain - rcoord = sqrt((x_cc(j)**2 + y_cc(k)**2 + z_cc(l)**2)) - theta = atan2(y_cc(k), x_cc(j)) - phi = atan2(sqrt(x_cc(j)**2 + y_cc(k)**2), z_cc(l)) - !spherical coord, assuming Rmax=1 - xi_sph = (rcoord**3 - R0ref**3 + 1_wp)**(1_wp/3_wp) - xi_cart(1) = xi_sph*sin(phi)*cos(theta) - xi_cart(2) = xi_sph*sin(phi)*sin(theta) - xi_cart(3) = xi_sph*cos(phi) - else - xi_cart(1) = x_cc(j) + if (hyperelasticity .and. .not. pre_stress) then + xi_cart(1) = x_cc(j) + if (p > 0) then xi_cart(2) = y_cc(k) xi_cart(3) = z_cc(l) + elseif (n > 0) then + xi_cart(2) = y_cc(k) end if ! assigning the reference map to the q_prim vector field diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 7f930b900e..670b4d067f 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -910,11 +910,6 @@ contains call MPI_INFO_CREATE(mpi_info_int, ierr) call MPI_INFO_SET(mpi_info_int, 'romio_ds_write', 'disable', ierr) - ! Option for UNIX file system (Hooke/Thomson) - ! WRITE(mpiiofs, '(A)') '/ufs_' - ! mpiiofs = TRIM(mpiiofs) - ! mpi_info_int = MPI_INFO_NULL - allocate (start_idx(1:num_dims)) #endif diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp index af627be391..32639b69f4 100644 --- a/src/pre_process/m_patches.fpp +++ b/src/pre_process/m_patches.fpp @@ -1523,7 +1523,8 @@ contains if (x_cc(i) - x_centroid >= 0 & .and. & - r - as(2)*Ps(2) - as(3)*Ps(3) - as(4)*Ps(4) - as(5)*Ps(5) - as(6)*Ps(6) - as(7)*Ps(7) - as(8)*Ps(8) - as(9)*Ps(9) <= radius .and. & + r - as(2)*Ps(2) - as(3)*Ps(3) - as(4)*Ps(4) - as(5)*Ps(5) - as(6)*Ps(6) & + - as(7)*Ps(7) - as(8)*Ps(8) - as(9)*Ps(9) <= radius .and. & patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & then call s_assign_patch_primitive_variables(patch_id, i, j, 0, & @@ -1531,7 +1532,8 @@ contains elseif (x_cc(i) - x_centroid < 0 & .and. & - r - as(2)*Ps(2) + as(3)*Ps(3) - as(4)*Ps(4) + as(5)*Ps(5) - as(6)*Ps(6) + as(7)*Ps(7) - as(8)*Ps(8) + as(9)*Ps(9) <= radius & + r - as(2)*Ps(2) + as(3)*Ps(3) - as(4)*Ps(4) + as(5)*Ps(5) - as(6)*Ps(6) & + + as(7)*Ps(7) - as(8)*Ps(8) + as(9)*Ps(9) <= radius & .and. & patch_icpp(patch_id)%alter_patch(patch_id_fp(i, j, 0))) & then diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 0503b40828..2a50584cfe 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -98,9 +98,9 @@ module m_cbc integer :: cbc_dir, cbc_loc !$acc declare create(dj, bcxb, bcxe, bcyb, bcye, bczb, bcze, cbc_dir, cbc_loc) - !! GRCBC inputs for subsonic inflow and outflow conditions consisting of - !! inflow velocities, pressure, density and void fraction as well as - !! outflow velocities and pressure + ! grCBC inputs for subsonic inflow and outflow conditions consisting of + ! inflow velocities, pressure, density and void fraction as well as + ! outflow velocities and pressure real(wp), allocatable, dimension(:) :: pres_in, pres_out, Del_in, Del_out real(wp), allocatable, dimension(:, :) :: vel_in, vel_out @@ -396,13 +396,13 @@ contains !$acc update device(bczb, bcze) end if - ! Allocate GRCBC inputs + ! Allocate grCBC inputs @:ALLOCATE(pres_in(1:num_dims), pres_out(1:num_dims)) @:ALLOCATE(Del_in(1:num_dims), Del_out(1:num_dims)) @:ALLOCATE(vel_in(1:num_dims, 1:num_dims), vel_out(1:num_dims, 1:num_dims)) @:ALLOCATE(alpha_rho_in(1:num_fluids, 1:num_dims), alpha_in(1:num_fluids, 1:num_dims)) - ! Assign and update GRCBC inputs + ! Assign and update grCBC inputs #:for CBC_DIR, XYZ in [(1, 'x'), (2, 'y'), (3, 'z')] if (${CBC_DIR}$ <= num_dims) then vel_in(${CBC_DIR}$, 1) = bc_${XYZ}$%vel_in(1) @@ -860,7 +860,7 @@ contains call s_compute_nonreflecting_subsonic_buffer_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) else if ((cbc_loc == -1 .and. bc${XYZ}$b == -7) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -7)) then call s_compute_nonreflecting_subsonic_inflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) - ! Add GRCBC for Subsonic Inflow + ! Add grCBC for Subsonic Inflow if (bc_${XYZ}$%grcbc_in) then !$acc loop seq do i = 2, momxb @@ -880,11 +880,11 @@ contains end if else if ((cbc_loc == -1 .and. bc${XYZ}$b == -8) .or. (cbc_loc == 1 .and. bc${XYZ}$e == -8)) then call s_compute_nonreflecting_subsonic_outflow_L(lambda, L, rho, c, mf, dalpha_rho_ds, dpres_ds, dvel_ds, dadv_ds) - ! Add GRCBC for Subsonic Outflow (Pressure) + ! Add grCBC for Subsonic Outflow (Pressure) if (bc_${XYZ}$%grcbc_out) then L(advxe) = c*(1._wp - Ma)*(pres - pres_out(${CBC_DIR}$))/Del_out(${CBC_DIR}$) - ! Add GRCBC for Subsonic Outflow (Normal Velocity) + ! Add grCBC for Subsonic Outflow (Normal Velocity) if (bc_${XYZ}$%grcbc_vel_out) then L(advxe) = L(advxe) + rho*c**2._wp*(1._wp - Ma)*(vel(dir_idx(1)) + vel_out(${CBC_DIR}$, dir_idx(1))*sign(1, cbc_loc))/Del_out(${CBC_DIR}$) end if @@ -1536,7 +1536,7 @@ contains ! Deallocating the cell-width distribution in the s-direction @:DEALLOCATE(ds) - ! Deallocating GRCBC inputs + ! Deallocating grCBC inputs @:DEALLOCATE(vel_in, vel_out, pres_in, pres_out, Del_in, Del_out, alpha_rho_in, alpha_in) ! Deallocating CBC Coefficients in x-direction diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index 6c819d2891..f18619723c 100644 --- a/src/simulation/m_checker.fpp +++ b/src/simulation/m_checker.fpp @@ -25,13 +25,12 @@ contains subroutine s_check_inputs call s_check_inputs_compilers - call s_check_inputs_weno call s_check_inputs_riemann_solver call s_check_inputs_time_stepping call s_check_inputs_model_eqns call s_check_inputs_acoustic_src - call s_check_inputs_hypoelasticity + call s_check_inputs_elasticity call s_check_inputs_bubbles_euler call s_check_inputs_bubbles_lagrange call s_check_inputs_adapt_dt @@ -123,7 +122,7 @@ contains @:PROHIBIT(model_eqns == 3 .and. wave_speeds /= 1, "6-equation model (model_eqns = 3) requires wave_speeds = 1") end subroutine s_check_inputs_model_eqns - !> Checks constraints for GRCBC + !> Checks constraints for grCBC subroutine s_check_inputs_grcbc #:for DIR in ['x', 'y', 'z'] @:PROHIBIT(bc_${DIR}$%grcbc_in .and. (bc_${DIR}$%beg /= -7 .and. bc_${DIR}$%end /= -7), "Subsonic Inflow requires bc = -7") @@ -259,9 +258,11 @@ contains end subroutine s_check_inputs_acoustic_src - !> Checks constraints on hypoelasticity parameters - subroutine s_check_inputs_hypoelasticity - @:PROHIBIT(hypoelasticity .and. riemann_solver /= 1, "hypoelasticity requires HLL Riemann solver (riemann_solver = 1)") + !> Checks constraints on elasticity parameters + subroutine s_check_inputs_elasticity + @:PROHIBIT(hyperelasticity .and. hyper_model == dflt_int) + @:PROHIBIT((hypoelasticity .or. hyperelasticity) .and. fd_order == dflt_int, & + "fd_order must be set for hypoelasticity or hyperelasticity") end subroutine !> Checks constraints on bubble parameters diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index eb9a1924ac..3b48bba713 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1084,7 +1084,7 @@ contains real(wp) :: nondim_time !< Non-dimensional time real(wp) :: tmp !< - !! Temporary variable to store quantity for mpi_allreduce + !! Temporary variable to store quantity for mpi_allreduce integer :: npts !< Number of included integral points real(wp) :: rad, thickness !< For integral quantities diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 26f435bc07..a5709d404c 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -710,7 +710,7 @@ contains integral(i)%ymax = dflt_real end do - ! GRCBC flags + ! grCBC flags #:for dir in {'x', 'y', 'z'} bc_${dir}$%grcbc_in = .false. bc_${dir}$%grcbc_out = .false. @@ -1086,6 +1086,27 @@ contains ! END: Volume Fraction Model + elasticity = hypoelasticity .or. hyperelasticity + + if (elasticity) then + ! creates stress indices for both hypo and hyperelasticity + stress_idx%beg = sys_size + 1 + stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2 + ! number of distinct stresses is 1 in 1D, 3 in 2D, 6 in 3D + sys_size = stress_idx%end + end if + + if (hyperelasticity) then + ! number of entries in the symmetric btensor plus the jacobian + b_size = (num_dims*(num_dims + 1))/2 + 1 + ! storing the jacobian in the last entry + tensor_size = num_dims**2 + 1 + xi_idx%beg = sys_size + 1 + xi_idx%end = sys_size + num_dims + ! adding three more equations for the \xi field and the elastic energy + sys_size = xi_idx%end + 1 + end if + if (chemistry) then species_idx%beg = sys_size + 1 species_idx%end = sys_size + num_species @@ -1256,11 +1277,6 @@ contains call MPI_INFO_CREATE(mpi_info_int, ierr) call MPI_INFO_SET(mpi_info_int, 'romio_ds_write', 'disable', ierr) - ! Option for UNIX file system (Hooke/Thomson) - ! WRITE(mpiiofs, '(A)') '/ufs_' - ! mpiiofs = TRIM(mpiiofs) - ! mpi_info_int = MPI_INFO_NULL - allocate (start_idx(1:num_dims)) #endif diff --git a/src/simulation/m_hyperelastic.fpp b/src/simulation/m_hyperelastic.fpp index 46addb7685..c005feeb78 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -116,53 +116,80 @@ contains alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l) alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l) end do + ! If in simulation, use acc mixture subroutines call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha_k, & alpha_rho_k, Re, j, k, l, G, Gs) rho = max(rho, sgm_eps) G = max(G, sgm_eps) - !if ( G <= verysmall ) G_K = 0_wp + !if ( G <= verysmall ) G_K = 0._wp if (G > verysmall) then !$acc loop seq - do i = 1, tensor_size - tensora(i) = 0_wp + do i = 1, tensor_size - 1 + tensora(i) = tensorb(i)/tensorb(tensor_size) end do + ! STEP 1: computing the grad_xi tensor using finite differences ! grad_xi definition / organization - ! number for the tensor 1-3: dxix_dx, dxiy_dx, dxiz_dx - ! 4-6 : dxix_dy, dxiy_dy, dxiz_dy - ! 7-9 : dxix_dz, dxiy_dz, dxiz_dz + ! tensora(1,2,3): dxix_dx, dxiy_dx, dxiz_dx + ! tensora(4,5,6): dxix_dy, dxiy_dy, dxiz_dy + ! tensora(7,8,9): dxix_dz, dxiy_dz, dxiz_dz !$acc loop seq do r = -fd_number, fd_number ! derivatives in the x-direction tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) - tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) - tensora(3) = tensora(3) + q_prim_vf(xiend)%sf(j + r, k, l)*fd_coeff_x(r, j) - ! derivatives in the y-direction - tensora(4) = tensora(4) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) - tensora(5) = tensora(5) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) - tensora(6) = tensora(6) + q_prim_vf(xiend)%sf(j, k + r, l)*fd_coeff_y(r, k) - ! derivatives in the z-direction - tensora(7) = tensora(7) + q_prim_vf(xibeg)%sf(j, k, l + r)*fd_coeff_z(r, l) - tensora(8) = tensora(8) + q_prim_vf(xibeg + 1)%sf(j, k, l + r)*fd_coeff_z(r, l) - tensora(9) = tensora(9) + q_prim_vf(xiend)%sf(j, k, l + r)*fd_coeff_z(r, l) + if (n > 0) then + ! derivatives in the x-direction + tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) + ! derivatives in the y-direction + tensora(3) = tensora(3) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(4) = tensora(4) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) + end if + if (p > 0) then + ! derivatives in the x-direction + tensora(2) = tensora(2) + q_prim_vf(xibeg + 1)%sf(j + r, k, l)*fd_coeff_x(r, j) + tensora(3) = tensora(3) + q_prim_vf(xiend)%sf(j + r, k, l)*fd_coeff_x(r, j) + ! derivatives in the y-direction + tensora(4) = tensora(4) + q_prim_vf(xibeg)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(5) = tensora(5) + q_prim_vf(xibeg + 1)%sf(j, k + r, l)*fd_coeff_y(r, k) + tensora(6) = tensora(6) + q_prim_vf(xiend)%sf(j, k + r, l)*fd_coeff_y(r, k) + ! derivatives in the z-direction + tensora(7) = tensora(7) + q_prim_vf(xibeg)%sf(j, k, l + r)*fd_coeff_z(r, l) + tensora(8) = tensora(8) + q_prim_vf(xibeg + 1)%sf(j, k, l + r)*fd_coeff_z(r, l) + tensora(9) = tensora(9) + q_prim_vf(xiend)%sf(j, k, l + r)*fd_coeff_z(r, l) + end if end do - ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse - tensorb(1) = tensora(5)*tensora(9) - tensora(6)*tensora(8) - tensorb(2) = -(tensora(2)*tensora(9) - tensora(3)*tensora(8)) - tensorb(3) = tensora(2)*tensora(6) - tensora(3)*tensora(5) - tensorb(4) = -(tensora(4)*tensora(9) - tensora(6)*tensora(7)) - tensorb(5) = tensora(1)*tensora(9) - tensora(3)*tensora(7) - tensorb(6) = -(tensora(1)*tensora(6) - tensora(4)*tensora(3)) - tensorb(7) = tensora(4)*tensora(8) - tensora(5)*tensora(7) - tensorb(8) = -(tensora(1)*tensora(8) - tensora(2)*tensora(7)) - tensorb(9) = tensora(1)*tensora(5) - tensora(2)*tensora(4) - - ! STEP 2b: computing the determinant of the grad_xi tensor - tensorb(tensor_size) = tensora(1)*(tensora(5)*tensora(9) - tensora(6)*tensora(8)) & - - tensora(2)*(tensora(4)*tensora(9) - tensora(6)*tensora(7)) & - + tensora(3)*(tensora(4)*tensora(8) - tensora(5)*tensora(7)) + + ! STEP 2a: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1) + ! STEP 2b: computing the inverse of the grad_xi tensor + tensorb(1) = 1._wp/(tensora(1)**2) + if (n > 0) then + ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse + tensorb(1) = tensora(4) + tensorb(2) = -tensora(3) + tensorb(3) = -tensora(2) + tensorb(4) = tensora(1) + ! STEP 2b: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1)*tensora(4) - tensora(2)*tensora(3) + end if + if (p > 0) then + ! STEP 2a: computing the adjoint of the grad_xi tensor for the inverse + tensorb(1) = tensora(5)*tensora(9) - tensora(6)*tensora(8) + tensorb(2) = -(tensora(2)*tensora(9) - tensora(3)*tensora(8)) + tensorb(3) = tensora(2)*tensora(6) - tensora(3)*tensora(5) + tensorb(4) = -(tensora(4)*tensora(9) - tensora(6)*tensora(7)) + tensorb(5) = tensora(1)*tensora(9) - tensora(3)*tensora(7) + tensorb(6) = -(tensora(1)*tensora(6) - tensora(4)*tensora(3)) + tensorb(7) = tensora(4)*tensora(8) - tensora(5)*tensora(7) + tensorb(8) = -(tensora(1)*tensora(8) - tensora(2)*tensora(7)) + tensorb(9) = tensora(1)*tensora(5) - tensora(2)*tensora(4) + ! STEP 2b: computing the determinant of the grad_xi tensor + tensorb(tensor_size) = tensora(1)*(tensora(5)*tensora(9) - tensora(6)*tensora(8)) & + - tensora(2)*(tensora(4)*tensora(9) - tensora(6)*tensora(7)) & + + tensora(3)*(tensora(4)*tensora(8) - tensora(5)*tensora(7)) + end if if (tensorb(tensor_size) > verysmall) then ! STEP 2c: computing the inverse of grad_xi tensor = F @@ -173,35 +200,53 @@ contains end do ! STEP 2d: computing the J = det(F) = 1/det(\grad{\xi}) - tensorb(tensor_size) = 1_wp/tensorb(tensor_size) - - ! STEP 3: computing F transpose F - tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensora(3)**2 - tensorb(5) = tensora(4)**2 + tensora(5)**2 + tensora(6)**2 - tensorb(9) = tensora(7)**2 + tensora(8)**2 + tensora(9)**2 - tensorb(2) = tensora(1)*tensora(4) + tensora(2)*tensora(5) + tensora(3)*tensora(6) - tensorb(3) = tensora(1)*tensora(7) + tensora(2)*tensora(8) + tensora(3)*tensora(9) - tensorb(6) = tensora(4)*tensora(7) + tensora(5)*tensora(8) + tensora(6)*tensora(9) - ! STEP 4: update the btensor, this is consistent with Riemann solvers - #:for BIJ, TXY in [(1,1),(2,2),(3,5),(4,3),(5,6),(6,9)] - btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$) - #:endfor - ! store the determinant at the last entry of the btensor + tensorb(tensor_size) = 1._wp/tensorb(tensor_size) + ! STEP 3: update the btensor, this is consistent with Riemann solvers + ! \b_xx + btensor%vf(1)%sf(j, k, l) = tensorb(1) + if (n > 0) then + ! STEP 2e: override adjoint (tensorb) to be F transpose F + tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensorb(4) = tensora(3)**2 + tensora(4)**2 + tensorb(2) = tensora(1)*tensora(3) + tensora(2)*tensora(4) + tensorb(3) = tensorb(2) !tensora(3)*tensora(1) + tensora(4)*tensora(2) + ! STEP 3: update the btensor, this is consistent with Riemann solvers + #:for BIJ, TXY in [(1,1),(2,2),(3,4)] + btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$) + #:endfor + end if + if (p > 0) then + ! STEP 2e: computing F transpose F + tensorb(1) = tensora(1)**2 + tensora(2)**2 + tensora(3)**2 + tensorb(5) = tensora(4)**2 + tensora(5)**2 + tensora(6)**2 + tensorb(9) = tensora(7)**2 + tensora(8)**2 + tensora(9)**2 + tensorb(2) = tensora(1)*tensora(4) + tensora(2)*tensora(5) + tensora(3)*tensora(6) + tensorb(3) = tensora(1)*tensora(7) + tensora(2)*tensora(8) + tensora(3)*tensora(9) + tensorb(6) = tensora(4)*tensora(7) + tensora(5)*tensora(8) + tensora(6)*tensora(9) + ! STEP 3a: update the btensor, this is consistent with Riemann solvers + #:for BIJ, TXY in [(1,1),(2,2),(3,5),(4,3),(5,6),(6,9)] + btensor%vf(${BIJ}$)%sf(j, k, l) = tensorb(${TXY}$) + #:endfor + end if + + !STEP 3b: store the determinant at the last entry of the btensor btensor%vf(b_size)%sf(j, k, l) = tensorb(tensor_size) - ! STEP 5a: updating the Cauchy stress primitive scalar field + + ! STEP 4a: updating the Cauchy stress primitive scalar field if (hyper_model == 1) then call s_neoHookean_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) elseif (hyper_model == 2) then call s_Mooney_Rivlin_cauchy_solver(btensor%vf, q_prim_vf, G, j, k, l) end if - ! STEP 5b: updating the pressure field + + ! STEP 4b: updating the pressure field q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) - & G*q_prim_vf(xiend + 1)%sf(j, k, l)/gamma - ! STEP 5c: updating the Cauchy stress conservative scalar field + + ! STEP 4c: updating the Cauchy stress conservative scalar field !$acc loop seq do i = 1, b_size - 1 - q_cons_vf(strxb + i - 1)%sf(j, k, l) = & - rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) + q_cons_vf(strxb + i - 1)%sf(j, k, l) = rho*q_prim_vf(strxb + i - 1)%sf(j, k, l) end do end if end if @@ -209,10 +254,12 @@ contains end do end do !$acc end parallel loop + end subroutine s_hyperelastic_rmt_stress_update - !> The following subroutine handles the calculation of the btensor. - !! The calculation of the btensor takes qprimvf. + !> The following subroutine handles the calculation of the btensor + !! with a neo-Hookean material model. + !! The calculation of the btensor takes qprimvf. !! @param q_prim_vf Primitive variables !! @param btensor is the output !! calculate the grad_xi, grad_xi is a nxn tensor @@ -228,15 +275,28 @@ contains real(wp) :: trace real(wp), parameter :: f13 = 1._wp/3._wp - integer :: i !< Generic loop iterators + integer :: i ! tensor is the symmetric tensor & calculate the trace of the tensor - trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l) - + trace = btensor(1)%sf(j, k, l) ! calculate the deviatoric of the tensor - #:for IJ in [1,3,6] - btensor(${IJ}$)%sf(j, k, l) = btensor(${IJ}$)%sf(j, k, l) - f13*trace - #:endfor + btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace + if (n > 0) then + ! tensor is the symmetric tensor & calculate the trace of the tensor + trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + ! calculate the deviatoric of the tensor + btensor(1)%sf(j, k, l) = btensor(1)%sf(j, k, l) - f13*trace + btensor(3)%sf(j, k, l) = btensor(3)%sf(j, k, l) - f13*trace + end if + if (p > 0) then + ! tensor is the symmetric tensor & calculate the trace of the tensor + trace = btensor(1)%sf(j, k, l) + btensor(3)%sf(j, k, l) + btensor(6)%sf(j, k, l) + ! calculate the deviatoric of the tensor + #:for IJ in [1,3,6] + btensor(${IJ}$)%sf(j, k, l) = btensor(${IJ}$)%sf(j, k, l) - f13*trace + #:endfor + end if + ! dividing by the jacobian for neo-Hookean model ! setting the tensor to the stresses for riemann solver !$acc loop seq @@ -250,7 +310,8 @@ contains end subroutine s_neoHookean_cauchy_solver - !> The following subroutine handles the calculation of the btensor. + !> The following subroutine handles the calculation of the btensor + !! with a Mooney-Rivlin material model. !! The calculation of the btensor takes qprimvf. !! @param q_prim_vf Primitive variables !! @param btensor is the output @@ -267,7 +328,7 @@ contains real(wp) :: trace real(wp), parameter :: f13 = 1._wp/3._wp - integer :: i !< Generic loop iterators + integer :: i !TODO Make this 1D and 2D capable ! tensor is the symmetric tensor & calculate the trace of the tensor diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp index 5f691de016..f0fb5f74e9 100644 --- a/src/simulation/m_hypoelastic.fpp +++ b/src/simulation/m_hypoelastic.fpp @@ -11,6 +11,7 @@ module m_hypoelastic use m_global_parameters !< Definitions of the global parameters use m_finite_differences use m_helper + use m_mpi_proxy !< Message passing interface (MPI) module proxy implicit none diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 5be17a762b..099336cd52 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -120,7 +120,7 @@ contains & 'wave_speeds', 'avg_state', 'precision', 'bc_x%beg', 'bc_x%end', & & 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'fd_order', & & 'num_probes', 'num_integrals', 'bubble_model', 'thermal', & - & 'R0_type', 'num_source', 'relax_model', 'num_ibs', 'n_start' ] + & 'R0_type', 'num_source', 'relax_model', 'hyper_model', 'num_ibs', 'n_start' ] call MPI_BCAST(${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) #:endfor diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index e2f8754add..3523ffd500 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -673,7 +673,10 @@ contains call nvtxEndRange call nvtxStartRange("RHS-ELASTIC") - if (hyperelasticity) call s_hyperelastic_rmt_stress_update(q_cons_qp%vf, q_prim_qp%vf) + if (hyperelasticity) then + call s_hyperelastic_rmt_stress_update(q_cons_qp%vf, q_prim_qp%vf) + call s_populate_variables_buffers(q_prim_qp%vf, pb, mv) + end if call nvtxEndRange if (cfl_dt) then @@ -712,7 +715,7 @@ contains call nvtxStartRange("RHS-WENO") if (.not. surface_tension) then - ! Reconstruct densitiess + ! Reconstruct densities iv%beg = 1; iv%end = sys_size call s_reconstruct_cell_boundary_values( & q_prim_qp%vf(1:sys_size), & diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 130a2aadae..ba7d92cb20 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -487,7 +487,6 @@ contains !$acc loop seq do i = 1, 2 Re_R(i) = dflt_real - if (Re_size(i) > 0) Re_R(i) = 0._wp !$acc loop seq @@ -545,8 +544,8 @@ contains call get_mixture_energy_mass(T_L, Ys_L, E_L) call get_mixture_energy_mass(T_R, Ys_R, E_R) - E_L = rho_L*E_L + 5e-1*rho_L*vel_L_rms - E_R = rho_R*E_R + 5e-1*rho_R*vel_R_rms + E_L = rho_L*E_L + 0.5_wp*rho_L*vel_L_rms + E_R = rho_R*E_R + 0.5_wp*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R elseif (mhd .and. relativity) then @@ -588,8 +587,8 @@ contains H_L = (E_L + pres_L - pres_mag%L)/rho_L H_R = (E_R + pres_R - pres_mag%R)/rho_R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) else - E_L = gamma_L*pres_L + pi_inf_L + 5e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5e-1*rho_R*vel_R_rms + qv_R + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R end if @@ -604,12 +603,13 @@ contains G_R = G_R + alpha_R(i)*Gs(i) end do + !$acc loop seq do i = 1, strxe - strxb + 1 tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) ! Elastic contribution to energy if G large enough !TODO take out if statement if stable without - if ((G_L > 1000) .and. (G_R > 1000)) then + if ((G_L > verysmall) .and. (G_R > verysmall)) then E_L = E_L + (tau_e_L(i)*tau_e_L(i))/(4._wp*G_L) E_R = E_R + (tau_e_R(i)*tau_e_R(i))/(4._wp*G_R) ! Double for shear stresses @@ -621,37 +621,31 @@ contains end do end if - ! elastic energy update - !if ( hyperelasticity ) then - ! G_L = 0._wp - ! G_R = 0._wp - ! - ! !$acc loop seq - ! do i = 1, num_fluids - ! G_L = G_L + alpha_L(i)*Gs(i) - ! G_R = G_R + alpha_R(i)*Gs(i) - ! end do - ! ! Elastic contribution to energy if G large enough - ! if ((G_L > 1e-3_wp) .and. (G_R > 1e-3_wp)) then - ! E_L = E_L + G_L*qL_prim_rs${XYZ}$_vf(j, k, l, xiend + 1) - ! E_R = E_R + G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, xiend + 1) - ! !$acc loop seq - ! do i = 1, b_size-1 - ! tau_e_L(i) = G_L*qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) - ! tau_e_R(i) = G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) - ! end do - ! !$acc loop seq - ! do i = 1, b_size-1 - ! tau_e_L(i) = 0_wp - ! tau_e_R(i) = 0_wp - ! end do - ! !$acc loop seq - ! do i = 1, num_dims - ! xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) - ! xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) - ! end do - ! end if - !end if + ! energy adjustments for hyperelastic energy + if (hyperelasticity) then + !$acc loop seq + do i = 1, num_dims + xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) + xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) + end do + G_L = 0._wp; G_R = 0._wp; + !$acc loop seq + do i = 1, num_fluids + ! Mixture left and right shear modulus + G_L = G_L + alpha_L(i)*Gs(i) + G_R = G_R + alpha_R(i)*Gs(i) + end do + ! Elastic contribution to energy if G large enough + if (G_L > verysmall .and. G_R > verysmall) then + E_L = E_L + G_L*qL_prim_rs${XYZ}$_vf(j, k, l, xiend + 1) + E_R = E_R + G_R*qR_prim_rs${XYZ}$_vf(j + 1, k, l, xiend + 1) + end if + !$acc loop seq + do i = 1, b_size - 1 + tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) + tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) + end do + end if @:compute_average_state() @@ -712,36 +706,37 @@ contains (s_R - vel_R(dir_idx(1)))) & /(rho_L*(s_L - vel_L(dir_idx(1))) - & rho_R*(s_R - vel_R(dir_idx(1)))) + elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(dir_idx(1)) - & + vel_R(dir_idx(1)))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(dir_idx(1)) - c_L*Ms_L s_R = vel_R(dir_idx(1)) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) - xi_M = (5e-1_wp + sign(5e-1_wp, s_L)) & - + (5e-1_wp - sign(5e-1_wp, s_L)) & - *(5e-1_wp + sign(5e-1_wp, s_R)) - xi_P = (5e-1_wp - sign(5e-1_wp, s_R)) & - + (5e-1_wp - sign(5e-1_wp, s_L)) & - *(5e-1_wp + sign(5e-1_wp, s_R)) + xi_M = (0.5_wp + sign(0.5_wp, s_L)) & + + (0.5_wp - sign(0.5_wp, s_L)) & + *(0.5_wp + sign(0.5_wp, s_R)) + xi_P = (0.5_wp - sign(0.5_wp, s_R)) & + + (0.5_wp - sign(0.5_wp, s_L)) & + *(0.5_wp + sign(0.5_wp, s_R)) ! Mass if (.not. relativity) then @@ -849,7 +844,7 @@ contains - rho_R*vel_R(dir_idx(i)))) & /(s_M - s_P) end do - else if (hypoelasticity) then + else if (elasticity) then !$acc loop seq do i = 1, num_vels flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = & @@ -903,7 +898,7 @@ contains - s_P*vel_L(dir_idx(1))*(E_L + pres_L - ptilde_L) & + s_M*s_P*(E_L - E_R)) & /(s_M - s_P) - else if (hypoelasticity) then + else if (elasticity) then !TODO: simplify this so it's not split into 3 if (num_dims == 1) then flux_rs${XYZ}$_vf(j, k, l, E_idx) = & @@ -944,9 +939,9 @@ contains /(s_M - s_P) end if - ! Elastic Stresses + ! elastic stresses flux. if (hypoelasticity) then - do i = 1, strxe - strxb + 1 !TODO: this indexing may be slow + do i = 1, strxe - strxb + 1 flux_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) = & (s_M*(rho_R*vel_R(dir_idx(1)) & *tau_e_R(i)) & @@ -958,7 +953,19 @@ contains end do end if - ! Advection + ! reference map flux. + if (hyperelasticity) then + do i = 1, num_dims + flux_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) = & + (s_M*rho_R*vel_R(dir_idx(1))*xi_field_R(i) & + - s_P*rho_L*vel_L(dir_idx(1))*xi_field_L(i) & + + s_M*s_P*(rho_L*xi_field_L(i) & + - rho_R*xi_field_R(i))) & + /(s_M - s_P) + end do + end if + + ! advection flux. !$acc loop seq do i = advxb, advxe flux_rs${XYZ}$_vf(j, k, l, i) = & @@ -971,18 +978,6 @@ contains /(s_M - s_P) end do - ! Xi field - !if ( hyperelasticity ) then - ! do i = 1, num_dims - ! flux_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) = & - ! (s_M*rho_R*vel_R(dir_idx(1))*xi_field_R(i) & - ! - s_P*rho_L*vel_L(dir_idx(1))*xi_field_L(i) & - ! + s_M*s_P*(rho_L*xi_field_L(i) & - ! - rho_R*xi_field_R(i))) & - ! /(s_M - s_P) - ! end do - !end if - ! Div(U)? !$acc loop seq do i = 1, num_vels @@ -1282,7 +1277,7 @@ contains !$acc private(vel_L, vel_R, vel_K_Star, Re_L, Re_R, rho_avg, h_avg, gamma_avg, & !$acc s_L, s_R, s_S, vel_avg_rms, alpha_L, alpha_R, Ys_L, Ys_R, Xs_L, Xs_R, & !$acc Gamma_iL, Gamma_iR, Cp_iL, Cp_iR, Yi_avg, Phi_avg, h_iL, h_iR, h_avg_2, & - !$acc tau_e_L, tau_e_R, G_L, G_R, flux_ene_e, xi_field_L, xi_field_R) + !$acc tau_e_L, tau_e_R, G_L, G_R, flux_ene_e, xi_field_L, xi_field_R, xi_L, xi_R) do l = is3%beg, is3%end do k = is2%beg, is2%end do j = is1%beg, is1%end @@ -1369,7 +1364,6 @@ contains Re_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, E_idx + Re_idx(i, q))/Res(i, q) & + Re_L(i) end do - Re_L(i) = 1._wp/max(Re_L(i), sgm_eps) end do @@ -1390,18 +1384,17 @@ contains end do end if - E_L = gamma_L*pres_L + pi_inf_L + 5e-1_wp*rho_L*vel_L_rms + qv_L - - E_R = gamma_R*pres_R + pi_inf_R + 5e-1_wp*rho_R*vel_R_rms + qv_R + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R - ! ENERGY ADJUSTMENTS FOR HYPOELASTIC ENERGY + ! energy adjustments for hypoelastic energy if (hypoelasticity) then !$acc loop seq do i = 1, strxe - strxb + 1 tau_e_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, strxb - 1 + i) tau_e_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, strxb - 1 + i) end do - G_L = 0_wp; G_R = 0_wp + G_L = 0._wp; G_R = 0._wp !$acc loop seq do i = 1, num_fluids G_L = G_L + alpha_L(i)*Gs(i) @@ -1422,14 +1415,14 @@ contains end do end if - ! ENERGY ADJUSTMENTS FOR HYPERELASTIC ENERGY + ! energy adjustments for hyperelastic energy if (hyperelasticity) then !$acc loop seq do i = 1, num_dims xi_field_L(i) = qL_prim_rs${XYZ}$_vf(j, k, l, xibeg - 1 + i) xi_field_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, xibeg - 1 + i) end do - G_L = 0_wp; G_R = 0_wp; + G_L = 0._wp; G_R = 0._wp; !$acc loop seq do i = 1, num_fluids ! Mixture left and right shear modulus @@ -1471,7 +1464,7 @@ contains end do end if - ! COMPUTING THE DIRECT WAVE SPEEDS + ! computing the direct wave speeds if (wave_speeds == 1) then if (elasticity) then s_L = min(vel_L(dir_idx(1)) - sqrt(c_L*c_L + & @@ -1493,25 +1486,25 @@ contains end if elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(dir_idx(1)) - & + vel_R(dir_idx(1)))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(dir_idx(1)) - c_L*Ms_L s_R = vel_R(dir_idx(1)) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if ! follows Einfeldt et al. @@ -1525,8 +1518,8 @@ contains ! goes with numerical star velocity in x/y/z directions ! xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5e-1_wp + sign(0.5_wp, s_S)) - xi_P = (5e-1_wp - sign(0.5_wp, s_S)) + xi_M = (0.5_wp + sign(0.5_wp, s_S)) + xi_P = (0.5_wp - sign(0.5_wp, s_S)) ! goes with the numerical velocity in x/y/z directions ! xi_P/M (pressure) = min/max(0. sgn(1,sL/sR)) @@ -1543,11 +1536,12 @@ contains rho_Star = xi_M*(rho_L*(xi_MP*xi_L + 1._wp - xi_MP)) + & xi_P*(rho_R*(xi_PP*xi_R + 1._wp - xi_PP)) - vel_K_Star = vel_L(idx1)*(1_wp - xi_MP) + xi_MP*vel_R(idx1) + & + vel_K_Star = vel_L(idx1)*(1._wp - xi_MP) + xi_MP*vel_R(idx1) + & xi_MP*xi_PP*(s_S - vel_R(idx1)) - ! COMPUTING FLUXES - ! MASS FLUX. + ! computing fluxes + + ! mass flux. !$acc loop seq do i = 1, contxe flux_rs${XYZ}$_vf(j, k, l, i) = & @@ -1555,30 +1549,30 @@ contains xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*(vel_R(idx1) + s_P*(xi_R - 1._wp)) end do - ! MOMENTUM FLUX. + ! momentum flux. ! f = \rho u u - \sigma, q = \rho u, q_star = \xi * \rho*(s_star, v, w) !$acc loop seq do i = 1, num_dims idxi = dir_idx(i) flux_rs${XYZ}$_vf(j, k, l, contxe + idxi) = rho_Star*vel_K_Star* & - (dir_flg(idxi)*vel_K_Star + (1_wp - dir_flg(idxi))*(xi_M*vel_L(idxi) + xi_P*vel_R(idxi))) + dir_flg(idxi)*p_Star + (dir_flg(idxi)*vel_K_Star + (1._wp - dir_flg(idxi))*(xi_M*vel_L(idxi) + xi_P*vel_R(idxi))) + dir_flg(idxi)*p_Star end do - ! ENERGY FLUX. + ! energy flux. ! f = u*(E-\sigma), q = E, q_star = \xi*E+(s-u)(\rho s_star - \sigma/(s-u)) flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_star + p_Star)*vel_K_Star - ! ELASTICITY. Elastic shear stress additions for the momentum and energy flux + ! elasticity. elastic shear stress additions for the momentum and energy flux if (elasticity) then - flux_ene_e = 0_wp; + flux_ene_e = 0._wp; !$acc loop seq do i = 1, num_dims idxi = dir_idx(i) - ! MOMENTUM ELASTIC FLUX. + ! momentum elastic flux. flux_rs${XYZ}$_vf(j, k, l, contxe + idxi) = & flux_rs${XYZ}$_vf(j, k, l, contxe + idxi) & - xi_M*tau_e_L(dir_idx_tau(i)) - xi_P*tau_e_R(dir_idx_tau(i)) - ! ENERGY ELASTIC FLUX. + ! energy elastic flux. flux_ene_e = flux_ene_e - & xi_M*(vel_L(idxi)*tau_e_L(dir_idx_tau(i)) + & s_M*(xi_L*((s_S - vel_L(i))*(tau_e_L(dir_idx_tau(i))/(s_L - vel_L(i)))))) - & @@ -1588,7 +1582,7 @@ contains flux_rs${XYZ}$_vf(j, k, l, E_idx) = flux_rs${XYZ}$_vf(j, k, l, E_idx) + flux_ene_e end if - ! VOLUME FRACTION FLUX. + ! volume fraction flux. !$acc loop seq do i = advxb, advxe flux_rs${XYZ}$_vf(j, k, l, i) = & @@ -1596,7 +1590,7 @@ contains xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*s_S end do - ! SOURCE TERM FOR VOLUME FRACTION ADVECTION FLUX. + ! source term for volume fraction advection flux. !$acc loop seq do i = 1, num_dims idxi = dir_idx(i) @@ -1605,14 +1599,14 @@ contains xi_P*(vel_R(idxi) + dir_flg(idxi)*(s_S*(xi_PP*(xi_R - 1) + 1) - vel_R(idxi))) end do - ! INTERNAL ENERGIES ADVECTION FLUX. + ! internal energies advection flux. ! K-th pressure and velocity in preparation for the internal energy flux !$acc loop seq do i = 1, num_fluids - p_K_Star = xi_M*(xi_MP*((pres_L + pi_infs(i)/(1_wp + gammas(i)))* & - xi_L**(1_wp/gammas(i) + 1_wp) - pi_infs(i)/(1_wp + gammas(i)) - pres_L) + pres_L) + & - xi_P*(xi_PP*((pres_R + pi_infs(i)/(1_wp + gammas(i)))* & - xi_R**(1_wp/gammas(i) + 1_wp) - pi_infs(i)/(1_wp + gammas(i)) - pres_R) + pres_R) + p_K_Star = xi_M*(xi_MP*((pres_L + pi_infs(i)/(1._wp + gammas(i)))* & + xi_L**(1._wp/gammas(i) + 1._wp) - pi_infs(i)/(1._wp + gammas(i)) - pres_L) + pres_L) + & + xi_P*(xi_PP*((pres_R + pi_infs(i)/(1._wp + gammas(i)))* & + xi_R**(1._wp/gammas(i) + 1._wp) - pi_infs(i)/(1._wp + gammas(i)) - pres_R) + pres_R) flux_rs${XYZ}$_vf(j, k, l, i + intxb - 1) = & ((xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, i + advxb - 1) + xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, i + advxb - 1))* & @@ -1623,7 +1617,7 @@ contains flux_src_rs${XYZ}$_vf(j, k, l, advxb) = vel_src_rs${XYZ}$_vf(j, k, l, idx1) - ! HYPOELASTIC STRESS EVOLUTION FLUX. + ! hypoelastic stress evolution flux. if (hypoelasticity) then !$acc loop seq do i = 1, strxe - strxb + 1 @@ -1633,7 +1627,7 @@ contains end do end if - ! REFERENCE MAP FLUX. + ! reference map flux. if (hyperelasticity) then !$acc loop seq do i = 1, num_dims @@ -1645,7 +1639,7 @@ contains end do end if - ! SURFACE TENSION FLUX. need to check + ! surface tension flux. need to check if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & (xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, c_idx) + & @@ -1670,7 +1664,7 @@ contains ! Geometrical source of the void fraction(s) is zero !$acc loop seq do i = advxb, advxe - flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0_wp + flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp end do end if #:endif @@ -1678,7 +1672,7 @@ contains if (grid_geometry == 3) then !$acc loop seq do i = 1, sys_size - flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0_wp + flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp end do flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb - 1 + dir_idx(1)) = & flux_gsrc_rs${XYZ}$_vf(j, k, l, momxb - 1 + dir_idx(1)) - p_Star @@ -1686,7 +1680,6 @@ contains flux_gsrc_rs${XYZ}$_vf(j, k, l, momxe) = flux_rs${XYZ}$_vf(j, k, l, momxb + 1) end if #:endif - end do end do end do @@ -1750,9 +1743,9 @@ contains qv_R = qv_R + alpha_rho_R(i)*qvs(i) end do - E_L = gamma_L*pres_L + pi_inf_L + 5e-1_wp*rho_L*vel_L_rms + qv_L + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5e-1_wp*rho_R*vel_R_rms + qv_R + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -1782,25 +1775,25 @@ contains /(rho_L*(s_L - vel_L(dir_idx(1))) - & rho_R*(s_R - vel_R(dir_idx(1)))) elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(dir_idx(1)) - & + vel_R(dir_idx(1)))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(dir_idx(1)) - c_L*Ms_L s_R = vel_R(dir_idx(1)) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if ! follows Einfeldt et al. @@ -1814,8 +1807,8 @@ contains ! goes with numerical velocity in x/y/z directions ! xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5e-1_wp + sign(5e-1_wp, s_S)) - xi_P = (5e-1_wp - sign(5e-1_wp, s_S)) + xi_M = (0.5_wp + sign(0.5_wp, s_S)) + xi_P = (0.5_wp - sign(0.5_wp, s_S)) !$acc loop seq do i = 1, contxe @@ -2046,7 +2039,6 @@ contains !$acc loop seq do i = 1, 2 Re_R(i) = dflt_real - if (Re_size(i) > 0) Re_R(i) = 0._wp !$acc loop seq @@ -2054,15 +2046,13 @@ contains Re_R(i) = (1._wp - qR_prim_rs${XYZ}$_vf(j + 1, k, l, E_idx + Re_idx(i, q)))/Res(i, q) & + Re_R(i) end do - Re_R(i) = 1._wp/max(Re_R(i), sgm_eps) end do end if end if - E_L = gamma_L*pres_L + pi_inf_L + 5e-1_wp*rho_L*vel_L_rms - - E_R = gamma_R*pres_R + pi_inf_R + 5e-1_wp*rho_R*vel_R_rms + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -2165,14 +2155,14 @@ contains if ((ptilde_L /= ptilde_L) .or. (ptilde_R /= ptilde_R)) then end if - rho_avg = 5e-1_wp*(rho_L + rho_R) - H_avg = 5e-1_wp*(H_L + H_R) - gamma_avg = 5e-1_wp*(gamma_L + gamma_R) + rho_avg = 0.5_wp*(rho_L + rho_R) + H_avg = 0.5_wp*(H_L + H_R) + gamma_avg = 0.5_wp*(gamma_L + gamma_R) vel_avg_rms = 0._wp !$acc loop seq do i = 1, num_dims - vel_avg_rms = vel_avg_rms + (5e-1_wp*(vel_L(i) + vel_R(i)))**2._wp + vel_avg_rms = vel_avg_rms + (0.5_wp*(vel_L(i) + vel_R(i)))**2._wp end do end if @@ -2210,25 +2200,25 @@ contains /(rho_L*(s_L - vel_L(dir_idx(1))) - & rho_R*(s_R - vel_R(dir_idx(1)))) elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(dir_idx(1)) - & - vel_R(dir_idx(1)))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(dir_idx(1)) - & + vel_R(dir_idx(1)))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(dir_idx(1)) - c_L*Ms_L s_R = vel_R(dir_idx(1)) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(dir_idx(1)) + vel_R(dir_idx(1))) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if ! follows Einfeldt et al. @@ -2242,8 +2232,8 @@ contains ! goes with numerical velocity in x/y/z directions ! xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5e-1_wp + sign(5e-1_wp, s_S)) - xi_P = (5e-1_wp - sign(5e-1_wp, s_S)) + xi_M = (0.5_wp + sign(0.5_wp, s_S)) + xi_P = (0.5_wp - sign(0.5_wp, s_S)) if (low_Mach == 1) then @:compute_low_Mach_correction() @@ -2570,20 +2560,19 @@ contains call get_mixture_energy_mass(T_L, Ys_L, E_L) call get_mixture_energy_mass(T_R, Ys_R, E_R) - E_L = rho_L*E_L + 5e-1*rho_L*vel_L_rms - E_R = rho_R*E_R + 5e-1*rho_R*vel_R_rms + E_L = rho_L*E_L + 0.5_wp*rho_L*vel_L_rms + E_R = rho_R*E_R + 0.5_wp*rho_R*vel_R_rms H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R else - E_L = gamma_L*pres_L + pi_inf_L + 5e-1*rho_L*vel_L_rms + qv_L - - E_R = gamma_R*pres_R + pi_inf_R + 5e-1*rho_R*vel_R_rms + qv_R + E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R end if - ! ENERGY ADJUSTMENTS FOR HYPOELASTIC ENERGY + ! energy adjustments for hypoelastic energy if (hypoelasticity) then !$acc loop seq do i = 1, strxe - strxb + 1 @@ -2612,7 +2601,7 @@ contains end do end if - ! ENERGY ADJUSTMENTS FOR HYPERELASTIC ENERGY + ! energy adjustments for hyperelastic energy if (hyperelasticity) then !$acc loop seq do i = 1, num_dims @@ -2687,25 +2676,25 @@ contains end if elseif (wave_speeds == 2) then - pres_SL = 5e-1_wp*(pres_L + pres_R + rho_avg*c_avg* & - (vel_L(idx1) - & - vel_R(idx1))) + pres_SL = 0.5_wp*(pres_L + pres_R + rho_avg*c_avg* & + (vel_L(idx1) - & + vel_R(idx1))) pres_SR = pres_SL - Ms_L = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_L)/(1._wp + gamma_L))* & + Ms_L = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_L)/(1._wp + gamma_L))* & (pres_SL/pres_L - 1._wp)*pres_L/ & ((pres_L + pi_inf_L/(1._wp + gamma_L))))) - Ms_R = max(1._wp, sqrt(1._wp + ((5e-1_wp + gamma_R)/(1._wp + gamma_R))* & + Ms_R = max(1._wp, sqrt(1._wp + ((0.5_wp + gamma_R)/(1._wp + gamma_R))* & (pres_SR/pres_R - 1._wp)*pres_R/ & ((pres_R + pi_inf_R/(1._wp + gamma_R))))) s_L = vel_L(idx1) - c_L*Ms_L s_R = vel_R(idx1) + c_R*Ms_R - s_S = 5e-1_wp*((vel_L(idx1) + vel_R(idx1)) + & - (pres_L - pres_R)/ & - (rho_avg*c_avg)) + s_S = 0.5_wp*((vel_L(idx1) + vel_R(idx1)) + & + (pres_L - pres_R)/ & + (rho_avg*c_avg)) end if ! follows Einfeldt et al. @@ -2719,11 +2708,12 @@ contains ! goes with numerical velocity in x/y/z directions ! xi_P/M = 0.5 +/m sgn(0.5,s_star) - xi_M = (5e-1_wp + sign(5e-1_wp, s_S)) - xi_P = (5e-1_wp - sign(5e-1_wp, s_S)) + xi_M = (0.5_wp + sign(0.5_wp, s_S)) + xi_P = (0.5_wp - sign(0.5_wp, s_S)) + + ! Computing the hllc fluxes - ! COMPUTING THE HLLC FLUXES - ! MASS FLUX. + ! mass flux. if (low_Mach == 1) then @:compute_low_Mach_correction() else @@ -2739,7 +2729,7 @@ contains *(vel_R(idx1) + s_P*(xi_R - 1._wp)) end do - ! MOMENTUM FLUX. + ! momentum flux. ! f = \rho u u - \sigma, q = \rho u, q_star = \xi * \rho*(s_star, v, w) !$acc loop seq do i = 1, num_dims @@ -2760,7 +2750,7 @@ contains + (s_M/s_L)*(s_P/s_R)*dir_flg(idxi)*pcorr end do - ! ENERGY FLUX. + ! energy flux. ! f = u*(E-\sigma), q = E, q_star = \xi*E+(s-u)(\rho s_star - \sigma/(s-u)) flux_rs${XYZ}$_vf(j, k, l, E_idx) = & xi_M*(vel_L(idx1)*(E_L + pres_L) + & @@ -2773,17 +2763,19 @@ contains (s_R - vel_R(idx1)))) - E_R)) & + (s_M/s_L)*(s_P/s_R)*pcorr*s_S - ! ELASTICITY. Elastic shear stress additions for the momentum and energy flux + ! elasticity. elastic shear stress additions for the momentum and energy flux if (elasticity) then flux_ene_e = 0_wp !$acc loop seq do i = 1, num_dims idxi = dir_idx(i) - ! MOMENTUM ELASTIC FLUX. + + ! momentum elastic flux flux_rs${XYZ}$_vf(j, k, l, contxe + idxi) = & flux_rs${XYZ}$_vf(j, k, l, contxe + idxi) & - xi_M*tau_e_L(dir_idx_tau(i)) - xi_P*tau_e_R(dir_idx_tau(i)) - ! ENERGY ELASTIC FLUX. + + ! energy elastic flux flux_ene_e = flux_ene_e - & xi_M*(vel_L(idxi)*tau_e_L(dir_idx_tau(i)) + & s_M*(xi_L*((s_S - vel_L(i))*(tau_e_L(dir_idx_tau(i))/(s_L - vel_L(i)))))) - & @@ -2793,7 +2785,7 @@ contains flux_rs${XYZ}$_vf(j, k, l, E_idx) = flux_rs${XYZ}$_vf(j, k, l, E_idx) + flux_ene_e end if - ! HYPOELASTIC STRESS EVOLUTION FLUX. + ! hypoelastic stress evolution flux if (hypoelasticity) then !$acc loop seq do i = 1, strxe - strxb + 1 @@ -2803,7 +2795,7 @@ contains end do end if - ! VOLUME FRACTION FLUX. + ! volume fraction flux !$acc loop seq do i = advxb, advxe flux_rs${XYZ}$_vf(j, k, l, i) = & @@ -2813,7 +2805,7 @@ contains *(vel_R(idx1) + s_P*(xi_R - 1._wp)) end do - ! VOLUME FRACTION SOURCE FLUX. + ! volume fraction source flux !$acc loop seq do i = 1, num_dims idxi = dir_idx(i) @@ -2826,7 +2818,7 @@ contains s_P*(xi_R - 1._wp)) end do - ! REFERENCE MAP FLUX. + ! reference map flux if (hyperelasticity) then !$acc loop seq do i = 1, num_dims @@ -3956,8 +3948,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j + 1, k, l)) tau_Re(1, 1) = (4._wp/3._wp)*dvel_avg_dx(1)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -3982,8 +3974,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dx(1)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -4010,18 +4002,18 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j + 1, k, l)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j + 1, k, l)) !$acc loop seq do i = 1, 2 dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j + 1, k, l)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j + 1, k, l)) end do - dvel_avg_dx(2) = 5e-1_wp*(dvelL_dx_vf(2)%sf(j, k, l) & - + dvelR_dx_vf(2)%sf(j + 1, k, l)) + dvel_avg_dx(2) = 0.5_wp*(dvelL_dx_vf(2)%sf(j, k, l) & + + dvelR_dx_vf(2)%sf(j + 1, k, l)) tau_Re(1, 1) = -(2._wp/3._wp)*(dvel_avg_dy(2) + & avg_vel(2)/y_cc(k))/ & @@ -4052,11 +4044,11 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j + 1, k, l)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j + 1, k, l)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j + 1, k, l)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j + 1, k, l)) tau_Re(1, 1) = (dvel_avg_dy(2) + & avg_vel(2)/y_cc(k))/ & @@ -4087,12 +4079,12 @@ contains !$acc loop seq do i = 1, 3, 2 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j + 1, k, l)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j + 1, k, l)) end do - dvel_avg_dx(3) = 5e-1_wp*(dvelL_dx_vf(3)%sf(j, k, l) & - + dvelR_dx_vf(3)%sf(j + 1, k, l)) + dvel_avg_dx(3) = 0.5_wp*(dvelL_dx_vf(3)%sf(j, k, l) & + + dvelR_dx_vf(3)%sf(j + 1, k, l)) tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dz(3)/y_cc(k)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -4125,8 +4117,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j + 1, k, l)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dz(3)/y_cc(k)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -4156,19 +4148,19 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k + 1, l)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j, k + 1, l)) !$acc loop seq do i = 1, 2 dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dx_vf(i)%sf(j, k, l) & + + dvelR_dx_vf(i)%sf(j, k + 1, l)) dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j, k + 1, l)) end do @@ -4205,14 +4197,14 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k + 1, l)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j, k + 1, l)) - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k + 1, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j, k + 1, l)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k + 1, l)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j, k + 1, l)) tau_Re(2, 2) = (dvel_avg_dx(1) + dvel_avg_dy(2) + & avg_vel(2)/y_cb(k))/ & @@ -4240,18 +4232,18 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(3) = 5e-1_wp*(velL_vf(3)%sf(j, k, l) & - + velR_vf(3)%sf(j, k + 1, l)) + avg_vel(3) = 0.5_wp*(velL_vf(3)%sf(j, k, l) & + + velR_vf(3)%sf(j, k + 1, l)) !$acc loop seq do i = 2, 3 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j, k + 1, l)) end do - dvel_avg_dy(3) = 5e-1_wp*(dvelL_dy_vf(3)%sf(j, k, l) & - + dvelR_dy_vf(3)%sf(j, k + 1, l)) + dvel_avg_dy(3) = 0.5_wp*(dvelL_dy_vf(3)%sf(j, k, l) & + + dvelR_dy_vf(3)%sf(j, k + 1, l)) tau_Re(2, 2) = -(2._wp/3._wp)*dvel_avg_dz(3)/y_cb(k)/ & Re_avg_rsy_vf(k, j, l, 1) @@ -4285,8 +4277,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k + 1, l)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j, k + 1, l)) tau_Re(2, 2) = dvel_avg_dz(3)/y_cb(k)/ & Re_avg_rsy_vf(k, j, l, 2) @@ -4317,28 +4309,28 @@ contains !$acc loop seq do i = 2, 3 - avg_vel(i) = 5e-1_wp*(velL_vf(i)%sf(j, k, l) & - + velR_vf(i)%sf(j, k, l + 1)) + avg_vel(i) = 0.5_wp*(velL_vf(i)%sf(j, k, l) & + + velR_vf(i)%sf(j, k, l + 1)) end do !$acc loop seq do i = 1, 3, 2 dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dx_vf(i)%sf(j, k, l) & + + dvelR_dx_vf(i)%sf(j, k, l + 1)) end do do i = 2, 3 dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j, k, l + 1)) end do !$acc loop seq do i = 1, 3 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j, k, l + 1)) end do tau_Re(3, 1) = (dvel_avg_dz(1)/y_cc(k) + dvel_avg_dx(3))/ & @@ -4380,17 +4372,17 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - avg_vel(2) = 5e-1_wp*(velL_vf(2)%sf(j, k, l) & - + velR_vf(2)%sf(j, k, l + 1)) + avg_vel(2) = 0.5_wp*(velL_vf(2)%sf(j, k, l) & + + velR_vf(2)%sf(j, k, l + 1)) - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k, l + 1)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j, k, l + 1)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k, l + 1)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j, k, l + 1)) - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k, l + 1)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j, k, l + 1)) tau_Re(3, 3) = (dvel_avg_dx(1) & + dvel_avg_dy(2) & @@ -4480,8 +4472,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j + 1, k, l)) tau_Re(1, 1) = (4._wp/3._wp)*dvel_avg_dx(1)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -4506,8 +4498,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j + 1, k, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dx(1)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -4537,12 +4529,12 @@ contains !$acc loop seq do i = 1, 2 dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j + 1, k, l)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j + 1, k, l)) end do - dvel_avg_dx(2) = 5e-1_wp*(dvelL_dx_vf(2)%sf(j, k, l) & - + dvelR_dx_vf(2)%sf(j + 1, k, l)) + dvel_avg_dx(2) = 0.5_wp*(dvelL_dx_vf(2)%sf(j, k, l) & + + dvelR_dx_vf(2)%sf(j + 1, k, l)) tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dy(2)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -4575,8 +4567,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j + 1, k, l)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dy(2)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -4606,12 +4598,12 @@ contains !$acc loop seq do i = 1, 3, 2 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j + 1, k, l)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j + 1, k, l)) end do - dvel_avg_dx(3) = 5e-1_wp*(dvelL_dx_vf(3)%sf(j, k, l) & - + dvelR_dx_vf(3)%sf(j + 1, k, l)) + dvel_avg_dx(3) = 0.5_wp*(dvelL_dx_vf(3)%sf(j, k, l) & + + dvelR_dx_vf(3)%sf(j + 1, k, l)) tau_Re(1, 1) = -(2._wp/3._wp)*dvel_avg_dz(3)/ & Re_avg_rsx_vf(j, k, l, 1) @@ -4643,8 +4635,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j + 1, k, l)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j + 1, k, l)) tau_Re(1, 1) = dvel_avg_dz(3)/ & Re_avg_rsx_vf(j, k, l, 2) @@ -4677,12 +4669,12 @@ contains do i = 1, 2 dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dx_vf(i)%sf(j, k, l) & + + dvelR_dx_vf(i)%sf(j, k + 1, l)) dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j, k + 1, l)) end do @@ -4718,11 +4710,11 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k + 1, l)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j, k + 1, l)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k + 1, l)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j, k + 1, l)) tau_Re(2, 2) = (dvel_avg_dx(1) + dvel_avg_dy(2))/ & Re_avg_rsy_vf(k, j, l, 2) @@ -4752,12 +4744,12 @@ contains !$acc loop seq do i = 2, 3 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k + 1, l)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j, k + 1, l)) end do - dvel_avg_dy(3) = 5e-1_wp*(dvelL_dy_vf(3)%sf(j, k, l) & - + dvelR_dy_vf(3)%sf(j, k + 1, l)) + dvel_avg_dy(3) = 0.5_wp*(dvelL_dy_vf(3)%sf(j, k, l) & + + dvelR_dy_vf(3)%sf(j, k + 1, l)) tau_Re(2, 2) = -(2._wp/3._wp)*dvel_avg_dz(3)/ & Re_avg_rsy_vf(k, j, l, 1) @@ -4790,8 +4782,8 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k + 1, l)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j, k + 1, l)) tau_Re(2, 2) = dvel_avg_dz(3)/ & Re_avg_rsy_vf(k, j, l, 2) @@ -4823,22 +4815,22 @@ contains !$acc loop seq do i = 1, 3, 2 dvel_avg_dx(i) = & - 5e-1_wp*(dvelL_dx_vf(i)%sf(j, k, l) & - + dvelR_dx_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dx_vf(i)%sf(j, k, l) & + + dvelR_dx_vf(i)%sf(j, k, l + 1)) end do !$acc loop seq do i = 2, 3 dvel_avg_dy(i) = & - 5e-1_wp*(dvelL_dy_vf(i)%sf(j, k, l) & - + dvelR_dy_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dy_vf(i)%sf(j, k, l) & + + dvelR_dy_vf(i)%sf(j, k, l + 1)) end do !$acc loop seq do i = 1, 3 dvel_avg_dz(i) = & - 5e-1_wp*(dvelL_dz_vf(i)%sf(j, k, l) & - + dvelR_dz_vf(i)%sf(j, k, l + 1)) + 0.5_wp*(dvelL_dz_vf(i)%sf(j, k, l) & + + dvelR_dz_vf(i)%sf(j, k, l + 1)) end do tau_Re(3, 1) = (dvel_avg_dz(1) + dvel_avg_dx(3))/ & @@ -4877,14 +4869,14 @@ contains do k = isy%beg, isy%end do j = isx%beg, isx%end - dvel_avg_dx(1) = 5e-1_wp*(dvelL_dx_vf(1)%sf(j, k, l) & - + dvelR_dx_vf(1)%sf(j, k, l + 1)) + dvel_avg_dx(1) = 0.5_wp*(dvelL_dx_vf(1)%sf(j, k, l) & + + dvelR_dx_vf(1)%sf(j, k, l + 1)) - dvel_avg_dy(2) = 5e-1_wp*(dvelL_dy_vf(2)%sf(j, k, l) & - + dvelR_dy_vf(2)%sf(j, k, l + 1)) + dvel_avg_dy(2) = 0.5_wp*(dvelL_dy_vf(2)%sf(j, k, l) & + + dvelR_dy_vf(2)%sf(j, k, l + 1)) - dvel_avg_dz(3) = 5e-1_wp*(dvelL_dz_vf(3)%sf(j, k, l) & - + dvelR_dz_vf(3)%sf(j, k, l + 1)) + dvel_avg_dz(3) = 0.5_wp*(dvelL_dz_vf(3)%sf(j, k, l) & + + dvelR_dz_vf(3)%sf(j, k, l + 1)) tau_Re(3, 3) = (dvel_avg_dx(1) & + dvel_avg_dy(2) & diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index 88d37f9674..263c8bf6fa 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -75,9 +75,9 @@ subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, pres = q_prim_vf(E_idx)%sf(j, k, l) - E = gamma*pres + pi_inf + 5e-1_wp*rho*vel_sum + qv + E = gamma*pres + pi_inf + 0.5_wp*rho*vel_sum + qv - ! ENERGY ADJUSTMENTS FOR HYPERELASTIC ENERGY + ! energy adjustments for hyperelastic energy if (hyperelasticity) then E = E + G*q_prim_vf(xiend + 1)%sf(j, k, l) end if @@ -246,7 +246,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) dy(k)/(abs(vel(2)) + c)) if (viscous) then - vcfl_dt = cfl_target*(min(dx(j), dy(k))**2._wp)/maxval((1/Re_l)/rho) + vcfl_dt = cfl_target*(min(dx(j), dy(k))**2._wp)/maxval((1._wp/Re_l)/rho) end if else @@ -254,7 +254,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) if (viscous) then - vcfl_dt = cfl_target*(dx(j)**2._wp)/minval(1/(rho*Re_l)) + vcfl_dt = cfl_target*(dx(j)**2._wp)/minval(1._wp/(rho*Re_l)) end if end if diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 9a5ed07fa1..3c7eb99c52 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -171,8 +171,9 @@ contains pi_fac, adv_n, adap_dt, bf_x, bf_y, bf_z, & k_x, k_y, k_z, w_x, w_y, w_z, p_x, p_y, p_z, & g_x, g_y, g_z, n_start, t_save, t_stop, & - cfl_adap_dt, cfl_const_dt, cfl_target, & - viscous, surface_tension, & + cfl_adap_dt, cfl_const_dt, cfl_target, & + viscous, surface_tension, & + hyperelasticity, hyper_model, R0ref, & bubbles_lagrange, lag_params, & rkck_adap_dt, rkck_tolerance, & hyperelasticity, R0ref, Bx0, powell diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 33836a9d9d..5522c28db8 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -852,12 +852,12 @@ contains call s_pressure_relaxation_procedure(q_cons_ts(1)%vf) end if + if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) + call nvtxStartRange("RHS-ELASTIC") if (hyperelasticity) call s_hyperelastic_rmt_stress_update(q_cons_ts(1)%vf, q_prim_vf) call nvtxEndRange - if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf) - if (ib) then if (qbmm .and. .not. polytropic) then call s_ibm_correct_state(q_cons_ts(1)%vf, q_prim_vf, pb_ts(1)%sf, mv_ts(1)%sf) @@ -1076,7 +1076,7 @@ contains start_rkck_step = .false. restart_rkck_step = .false. - ! FIRST TIME-STAGE + ! first time-stage RKstep = 1 rkck_time_tmp = mytime + rkck_c1*dt !$acc update device (rkck_time_tmp) @@ -1090,7 +1090,7 @@ contains if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step) if (restart_rkck_step) cycle - ! SECOND TIME-STAGE + ! second time-stage RKstep = 2 rkck_time_tmp = mytime + rkck_c2*dt !$acc update device (rkck_time_tmp) @@ -1104,7 +1104,7 @@ contains if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step) if (restart_rkck_step) cycle - ! THIRD TIME-STAGE + ! third time-stage RKstep = 3 rkck_time_tmp = mytime + rkck_c3*dt !$acc update device (rkck_time_tmp) @@ -1118,7 +1118,7 @@ contains if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step) if (restart_rkck_step) cycle - ! FOURTH TIME-STAGE + ! fourth time-stage RKstep = 4 rkck_time_tmp = mytime + rkck_c4*dt !$acc update device (rkck_time_tmp) @@ -1132,7 +1132,7 @@ contains if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step) if (restart_rkck_step) cycle - ! FIFTH TIME-STAGE + ! fifth time-stage RKstep = 5 rkck_time_tmp = mytime + rkck_c5*dt !$acc update device (rkck_time_tmp) @@ -1146,7 +1146,7 @@ contains if (lag_largestep > 0._wp) call s_compute_rkck_dt(lag_largestep, restart_rkck_step) if (restart_rkck_step) cycle - ! SIXTH TIME-STAGE + ! sixth time-stage RKstep = 6 rkck_time_tmp = mytime + rkck_c6*dt !$acc update device (rkck_time_tmp) @@ -1163,7 +1163,7 @@ contains dt_did = dt if (rkck_adap_dt) then - ! TRUNCATION ERROR + ! truncation error #ifdef DEBUG if (proc_rank == 0) print *, 'Computing truncation error (4th/5th RKCK)' #endif diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 65f0c5151b..3510827a64 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -253,6 +253,7 @@ def analytic(self): 'low_Mach': ParamType.INT, 'surface_tension': ParamType.LOG, 'viscous': ParamType.LOG, + 'hyper_model': ParamType.INT, 'bubbles_lagrange': ParamType.LOG, 'rkck_tolerance': ParamType.REAL, 'powell': ParamType.LOG, @@ -372,6 +373,7 @@ def analytic(self): 'flux_wrt': ParamType.LOG, 'E_wrt': ParamType.LOG, 'pres_wrt': ParamType.LOG, + 'tau_wrt': ParamType.LOG, 'alpha_wrt': ParamType.LOG, 'kappa_wrt': ParamType.LOG, 'gamma_wrt': ParamType.LOG, @@ -393,6 +395,7 @@ def analytic(self): 't_stop': ParamType.REAL, 'n_start': ParamType.INT, 'surface_tension': ParamType.LOG, + 'kymograph': ParamType.LOG, 'output_partial_domain': ParamType.LOG, 'bubbles_lagrange': ParamType.LOG, }) diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index c7d366af06..ed95837eb9 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -535,7 +535,7 @@ def alter_hypoelasticity(dimInfo): # Hypoelasticity checks for num_fluids in [1,2]: stack.push(f"Hypoelasticity -> {num_fluids} Fluid(s)", { - "hypoelasticity": 'T', "num_fluids": num_fluids, + 'hypoelasticity': 'T', 'num_fluids': num_fluids, 'riemann_solver': 1, 'fd_order': 4, 'fluid_pp(1)%gamma': 0.3, 'fluid_pp(1)%pi_inf': 7.8E+05, @@ -900,8 +900,8 @@ def foreach_example(): if path == "scaling": continue - # # List of currently broken examples -> currently attempting to fix! - brokenCases = ["2D_ibm_cfl_dt", "1D_sodHypo", "2D_viscous", "2D_laplace_pressure_jump", "2D_bubbly_steady_shock", "2D_advection", "2D_hardcodied_ic", "2D_ibm_multiphase", "2D_acoustic_broadband", "1D_inert_shocktube", "1D_reactive_shocktube", "2D_ibm_steady_shock", "3D_performance_test", "3D_ibm_stl_ellipsoid", "3D_sphbubcollapse", "2D_ibm_stl_wedge", "3D_ibm_stl_pyramid", "3D_ibm_bowshock", "3D_turb_mixing", "2D_mixing_artificial_Ma", "3D_lagrange_bubblescreen", "2D_triple_point"] + # List of currently broken examples -> currently attempting to fix! + brokenCases = ["2D_ibm_cfl_dt", "1D_sodHypo", "2D_viscous", "2D_laplace_pressure_jump", "2D_bubbly_steady_shock", "2D_advection", "2D_hardcodied_ic", "2D_ibm_multiphase", "2D_acoustic_broadband", "1D_inert_shocktube", "1D_reactive_shocktube", "2D_ibm_steady_shock", "3D_performance_test", "3D_ibm_stl_ellipsoid", "3D_sphbubcollapse", "2D_ibm_stl_wedge", "3D_ibm_stl_pyramid", "3D_ibm_bowshock", "3D_turb_mixing", "2D_mixing_artificial_Ma", "3D_lagrange_bubblescreen", "3D_hyper_bubingel", "3D_hyper_bubinwater", "2D_triple_point", "1D_hyper_impact_strong", "1D_hyper_impact_weak", "1D_hypo_impact_strong", "1D_hypo_impact_weak" ] if path in brokenCases: continue name = f"{path.split('_')[0]} -> Example -> {'_'.join(path.split('_')[1:])}" diff --git a/toolchain/templates/oscar.mako b/toolchain/templates/oscar.mako index 158a217cc1..f837553859 100644 --- a/toolchain/templates/oscar.mako +++ b/toolchain/templates/oscar.mako @@ -15,9 +15,8 @@ #SBATCH --account="${account}" % endif % if gpu: -#SBATCH --gpus-per-node=${tasks_per_node} -#SBATCH --mem=64G -#SBATCH --gpu-bind=closest +#SBATCH --gpu-bind=verbose,closest +#SBATCH --gres=gpu:v100-16:${tasks_per_node} % endif #SBATCH --output="${name}.out" #SBATCH --error="${name}.err" @@ -31,7 +30,7 @@ ${helpers.template_prologue()} ok ":) Loading modules:\n" -cd "${MFC_ROOTDIR}" +cd "${MFC_ROOT_DIR}" . ./mfc.sh load -c o -m ${'g' if gpu else 'c'} cd - > /dev/null echo @@ -42,9 +41,8 @@ echo % if not mpi: (set -x; ${profiler} "${target.get_install_binpath(case)}") % else: - (set -x; ${profiler} \ - mpirun -np ${nodes*tasks_per_node} \ - ${' '.join([f"'{x}'" for x in ARG('--') ])} \ + (set -x; ${profiler} \ + mpirun -np ${nodes*tasks_per_node} \ "${target.get_install_binpath(case)}") % endif