diff --git a/.github/workflows/lint-source.yml b/.github/workflows/lint-source.yml index cc6757697a..7b78d15254 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 5b3f54b202..c74f421a96 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -508,6 +508,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_1material/case.py b/examples/1D_hyper_1material/case.py new file mode 100644 index 0000000000..823021071f --- /dev/null +++ b/examples/1D_hyper_1material/case.py @@ -0,0 +1,107 @@ +#!/usr/bin/env python3 + +# 1D -> Hyperelasticity -> 1 Fluid(s) + +import math, json + +print(json.dumps({ + 'run_time_info': 'T', + 'm': 299, + 'n': 0, + 'p': 0, + 'dt': 1e-07, + 't_step_start': 0, + 't_step_stop': 50, + 't_step_save': 25, + 'num_patches': 3, + 'model_eqns': 3, + 'alt_soundspeed': 'F', + 'num_fluids': 1, + 'mpp_lim': 'F', + 'mixture_err': 'F', + 'time_stepper': 3, + 'weno_order': 5, + 'weno_eps': 1e-16, + 'mapped_weno': 'F', + 'null_weights': 'F', + 'mp_weno': 'F', + 'riemann_solver': 2, + 'wave_speeds': 1, + 'avg_state': 2, + 'format': 1, + 'precision': 2, + 'prim_vars_wrt': 'F', + 'parallel_io': 'F', + 'patch_icpp(1)%pres': 1000000.0, + 'patch_icpp(1)%alpha_rho(1)': 1000.0, + 'patch_icpp(1)%alpha(1)': 1.0, + 'patch_icpp(2)%pres': 100000.0, + 'patch_icpp(2)%alpha_rho(1)': 1000.0, + 'patch_icpp(2)%alpha(1)': 1.0, + 'patch_icpp(3)%pres': 500000.0, + 'patch_icpp(3)%alpha_rho(1)': 1000.0, + 'patch_icpp(3)%alpha(1)': 1.0, + 'fluid_pp(1)%gamma': 0.3, + 'fluid_pp(1)%pi_inf': 780000.0, + 'Ca': 0.9769178386380458, + 'Web': 13.927835051546392, + 'Re_inv': 0.009954269975623245, + 'pref': 101325.0, + 'rhoref': 1000.0, + 'bubble_model': 3, + 'polytropic': 'T', + 'polydisperse': 'F', + 'thermal': 3, + 'R0ref': 1e-05, + 'patch_icpp(1)%r0': 1, + 'patch_icpp(1)%v0': 0, + 'patch_icpp(2)%r0': 1, + 'patch_icpp(2)%v0': 0, + 'patch_icpp(3)%r0': 1, + 'patch_icpp(3)%v0': 0, + 'qbmm': 'F', + 'dist_type': 2, + 'poly_sigma': 0.3, + 'R0_type': 1, + 'sigR': 0.1, + 'sigV': 0.1, + 'rhoRV': 0.0, + 'x_domain%beg': 0.0, + 'x_domain%end': 1.0, + 'bc_x%beg': -3, + 'bc_x%end': -3, + 'patch_icpp(1)%geometry': 1, + 'patch_icpp(1)%x_centroid': 0.05, + 'patch_icpp(1)%length_x': 0.1, + 'patch_icpp(2)%x_centroid': 0.45, + 'patch_icpp(2)%length_x': 0.7, + 'patch_icpp(3)%x_centroid': 0.9, + 'patch_icpp(3)%length_x': 0.2, + 'patch_icpp(1)%vel(1)': 0.0, + 'patch_icpp(2)%geometry': 1, + 'patch_icpp(2)%vel(1)': 0.0, + 'patch_icpp(3)%geometry': 1, + 'patch_icpp(3)%vel(1)': 0.0, + 'hyperelasticity': 'T', + 'hyper_model': 1, + 'fd_order': 4, + 'patch_icpp(1)%tau_e(1)': 0.0, + 'patch_icpp(2)%tau_e(1)': 0.0, + 'patch_icpp(3)%tau_e(1)': 0.0, + 'fluid_pp(1)%G': 1000000000.0, + 'parallel_io' : 'T', + 'cons_vars_wrt' : 'T', + 'prim_vars_wrt': 'T', + 'alpha_rho_wrt(1)': 'T', + 'rho_wrt' : 'T', + 'mom_wrt(1)' : 'T', + 'vel_wrt(1)' : 'T', + 'E_wrt' : 'T', + 'pres_wrt' : 'T', + 'alpha_wrt(1)' : 'T', + 'gamma_wrt' : 'T', + 'heat_ratio_wrt' : 'T', + 'pi_inf_wrt' : 'T', + 'pres_inf_wrt' : 'T', + 'c_wrt' : 'T', +})) diff --git a/examples/1D_hyper_2materials/case.py b/examples/1D_hyper_2materials/case.py new file mode 100644 index 0000000000..4621386498 --- /dev/null +++ b/examples/1D_hyper_2materials/case.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python3 + +# 1D -> Hyperelasticity -> 2 Fluid(s) + +import math, json +import argparse + +print(json.dumps({ + 'run_time_info': 'T', + 'm': 299, + 'n': 0, + 'p': 0, + 'dt': 1e-07, + 't_step_start': 0, + 't_step_stop': 50, + 't_step_save': 25, + 'num_patches': 3, + 'model_eqns': 3, + 'alt_soundspeed': 'F', + 'num_fluids': 2, + 'mpp_lim': 'F', + 'mixture_err': 'F', + 'time_stepper': 3, + 'weno_order': 5, + 'weno_eps': 1e-16, + 'mapped_weno': 'F', + 'null_weights': 'F', + 'mp_weno': 'F', + 'riemann_solver': 2, + 'wave_speeds': 1, + 'avg_state': 2, + 'format': 1, + 'precision': 2, + 'prim_vars_wrt': 'F', + 'parallel_io': 'F', + 'patch_icpp(1)%pres': 1000000.0, + 'patch_icpp(1)%alpha_rho(1)': 900.0, + 'patch_icpp(1)%alpha(1)': 0.9, + 'patch_icpp(2)%pres': 100000.0, + 'patch_icpp(2)%alpha_rho(1)': 100, + 'patch_icpp(2)%alpha(1)': 0.1, + 'patch_icpp(3)%pres': 500000.0, + 'patch_icpp(3)%alpha_rho(1)': 900, + 'patch_icpp(3)%alpha(1)': 0.9, + 'fluid_pp(1)%gamma': 0.3, + 'fluid_pp(1)%pi_inf': 780000.0, +# 'Ca': 0.9769178386380458, +# 'Web': 13.927835051546392, +# 'Re_inv': 0.009954269975623245, + 'pref': 101325.0, + 'rhoref': 1000.0, +# 'bubble_model': 3, +# 'polytropic': 'T', +# 'polydisperse': 'F', +# 'thermal': 3, + 'R0ref': 1e-05, + 'patch_icpp(1)%r0': 1, + 'patch_icpp(1)%v0': 0, + 'patch_icpp(2)%r0': 1, + 'patch_icpp(2)%v0': 0, + 'patch_icpp(3)%r0': 1, + 'patch_icpp(3)%v0': 0, +# 'qbmm': 'F', +# 'dist_type': 2, +# 'poly_sigma': 0.3, +# 'R0_type': 1, +# 'sigR': 0.1, +# 'sigV': 0.1, +# 'rhoRV': 0.0, + 'x_domain%beg': 0.0, + 'x_domain%end': 1.0, + 'bc_x%beg': -3, + 'bc_x%end': -3, + 'patch_icpp(1)%geometry': 1, + 'patch_icpp(1)%x_centroid': 0.05, + 'patch_icpp(1)%length_x': 0.1, + 'patch_icpp(2)%x_centroid': 0.45, + 'patch_icpp(2)%length_x': 0.7, + 'patch_icpp(3)%x_centroid': 0.9, + 'patch_icpp(3)%length_x': 0.2, + 'patch_icpp(1)%vel(1)': 0.0, + 'patch_icpp(2)%geometry': 1, + 'patch_icpp(2)%vel(1)': 0.0, + 'patch_icpp(3)%geometry': 1, + 'patch_icpp(3)%vel(1)': 0.0, + 'hyperelasticity': 'T', + 'hyper_model': 1, + 'fd_order': 4, + 'patch_icpp(1)%tau_e(1)': 0.0, + 'patch_icpp(2)%tau_e(1)': 0.0, + 'patch_icpp(3)%tau_e(1)': 0.0, + 'fluid_pp(1)%G': 1000000000.0, + 'fluid_pp(2)%gamma': 0.3, + 'fluid_pp(2)%pi_inf': 780000.0, + 'patch_icpp(1)%alpha_rho(2)': 100, + 'patch_icpp(1)%alpha(2)': 0.1, + 'patch_icpp(2)%alpha_rho(2)': 900, + 'patch_icpp(2)%alpha(2)': 0.9, + 'patch_icpp(3)%alpha_rho(2)': 100, + 'patch_icpp(3)%alpha(2)': 0.1, + 'fluid_pp(2)%G': 0.0, #50000.0 + 'parallel_io' : 'T', 'cons_vars_wrt' : 'T', + 'prim_vars_wrt': 'T', 'alpha_rho_wrt(1)': 'T', + 'rho_wrt' : 'T', 'mom_wrt(1)' : 'T', + 'vel_wrt(1)' : 'T', 'E_wrt' : 'T', + 'pres_wrt' : 'T', 'alpha_wrt(1)' : 'T', + 'gamma_wrt' : 'T', 'heat_ratio_wrt' : 'T', + 'pi_inf_wrt' : 'T', 'pres_inf_wrt' : 'T', + 'c_wrt' : 'T', + })) 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/2D_hyper_1material/case.py b/examples/2D_hyper_1material/case.py new file mode 100644 index 0000000000..94890e1335 --- /dev/null +++ b/examples/2D_hyper_1material/case.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python3 + +# 2D -> Hyperelasticity -> 1 Fluid(s) + +import math, json + +print(json.dumps({ + 'run_time_info': 'T', + 'm': 49, + 'n': 39, + 'p': 0, + 'dt': 1e-08, #1e-06, +# 'cfl_adap_dt': 'T', +# 'cfl_target': 0.1, +# 'n_start': 0, +# 't_stop': 1e-06, +# 't_save': 1e-07, + 't_step_start': 0, + 't_step_stop': 50, + 't_step_save': 1, + 'num_patches': 3, + 'model_eqns': 3, + 'alt_soundspeed': 'F', + 'num_fluids': 1, + 'mpp_lim': 'F', + 'mixture_err': 'F', + 'time_stepper': 3, + 'weno_order': 5, + 'weno_eps': 1e-16, + 'mapped_weno': 'F', + 'null_weights': 'F', + 'mp_weno': 'F', + 'riemann_solver': 2, + 'wave_speeds': 1, + 'avg_state': 2, + 'format': 1, + 'precision': 2, + 'prim_vars_wrt': 'F', + 'parallel_io': 'F', + 'patch_icpp(1)%pres': 1000000.0, + 'patch_icpp(1)%alpha_rho(1)': 1000.0, + 'patch_icpp(1)%alpha(1)': 1.0, + 'patch_icpp(2)%pres': 100000.0, + 'patch_icpp(2)%alpha_rho(1)': 1000.0, + 'patch_icpp(2)%alpha(1)': 1.0, + 'patch_icpp(3)%pres': 500000.0, + 'patch_icpp(3)%alpha_rho(1)': 1000.0, + 'patch_icpp(3)%alpha(1)': 1.0, + 'fluid_pp(1)%gamma': 0.3, + 'fluid_pp(1)%pi_inf': 780000.0, +# 'Ca': 0.9769178386380458, +# 'Web': 13.927835051546392, +# 'Re_inv': 0.009954269975623245, + 'pref': 101325.0, + 'rhoref': 1000.0, +# 'bubble_model': 3, +# 'polytropic': 'T', +# 'polydisperse': 'F', +# 'thermal': 3, + 'R0ref': 1e-05, + 'patch_icpp(1)%r0': 1, + 'patch_icpp(1)%v0': 0, + 'patch_icpp(2)%r0': 1, + 'patch_icpp(2)%v0': 0, + 'patch_icpp(3)%r0': 1, + 'patch_icpp(3)%v0': 0, +# 'qbmm': 'F', +# 'dist_type': 2, +# 'poly_sigma': 0.3, +# 'R0_type': 1, +# 'sigR': 0.1, +# 'sigV': 0.1, +# 'rhoRV': 0.0, + 'x_domain%beg': 0.0, + 'x_domain%end': 1.0, + 'y_domain%beg': 0.0, + 'y_domain%end': 1.0, + 'bc_x%beg': -3, + 'bc_x%end': -3, + 'bc_y%beg': -3, + 'bc_y%end': -3, + 'patch_icpp(1)%geometry': 3, + 'patch_icpp(1)%y_centroid': 0.05, + 'patch_icpp(1)%length_y': 0.1, + 'patch_icpp(2)%y_centroid': 0.45, + 'patch_icpp(2)%length_y': 0.7, + 'patch_icpp(3)%y_centroid': 0.9, + 'patch_icpp(3)%length_y': 0.2, + 'patch_icpp(1)%x_centroid': 0.5, + 'patch_icpp(1)%length_x': 1, + 'patch_icpp(1)%vel(1)': 0.0, + 'patch_icpp(1)%vel(2)': 0.0, + 'patch_icpp(2)%geometry': 3, + 'patch_icpp(2)%x_centroid': 0.5, + 'patch_icpp(2)%length_x': 1, + 'patch_icpp(2)%vel(1)': 0.0, + 'patch_icpp(2)%vel(2)': 0.0, + 'patch_icpp(3)%geometry': 3, + 'patch_icpp(3)%x_centroid': 0.5, + 'patch_icpp(3)%length_x': 1, + 'patch_icpp(3)%vel(1)': 0.0, + 'patch_icpp(3)%vel(2)': 0.0, + 'hyperelasticity': 'T', + 'hyper_model': 1, + 'fd_order': 4, + 'patch_icpp(1)%tau_e(1)': 0.0, + 'patch_icpp(2)%tau_e(1)': 0.0, + 'patch_icpp(3)%tau_e(1)': 0.0, + 'fluid_pp(1)%G': 1000000000.0, + 'patch_icpp(1)%tau_e(2)': 0.0, + 'patch_icpp(1)%tau_e(3)': 0.0, + 'patch_icpp(2)%tau_e(2)': 0.0, + 'patch_icpp(2)%tau_e(3)': 0.0, + 'patch_icpp(3)%tau_e(2)': 0.0, + 'patch_icpp(3)%tau_e(3)': 0.0, + 'parallel_io' : 'T', 'cons_vars_wrt' : 'T', + 'prim_vars_wrt': 'T', 'alpha_rho_wrt(1)': 'T', + 'rho_wrt' : 'T', 'mom_wrt(1)' : 'T', + 'vel_wrt(1)' : 'T', 'E_wrt' : 'T', + 'pres_wrt' : 'T', 'alpha_wrt(1)' : 'T', + 'gamma_wrt' : 'T', 'heat_ratio_wrt' : 'T', + 'pi_inf_wrt' : 'T', 'pres_inf_wrt' : 'T', + 'c_wrt' : 'T', + })) + diff --git a/examples/2D_hyper_2materials/case.py b/examples/2D_hyper_2materials/case.py new file mode 100644 index 0000000000..0bc60ac510 --- /dev/null +++ b/examples/2D_hyper_2materials/case.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 + +# 2D -> Hyperelasticity -> 2 Fluid(s) + +import math, json + +print(json.dumps({ + 'run_time_info': 'T', + 'm': 49, + 'n': 39, + 'p': 0, + 'dt': 1e-06, + 't_step_start': 0, + 't_step_stop': 50, + 't_step_save': 25, + 'num_patches': 3, + 'model_eqns': 3, + 'alt_soundspeed': 'F', + 'num_fluids': 2, + 'mpp_lim': 'F', + 'mixture_err': 'F', + 'time_stepper': 3, + 'weno_order': 5, + 'weno_eps': 1e-16, + 'mapped_weno': 'F', + 'null_weights': 'F', + 'mp_weno': 'F', + 'riemann_solver': 2, + 'wave_speeds': 1, + 'avg_state': 2, + 'format': 1, + 'precision': 2, + 'prim_vars_wrt': 'F', + 'parallel_io': 'F', + 'patch_icpp(1)%pres': 1000000.0, + 'patch_icpp(1)%alpha_rho(1)': 900.0, + 'patch_icpp(1)%alpha(1)': 0.9, + 'patch_icpp(2)%pres': 100000.0, + 'patch_icpp(2)%alpha_rho(1)': 100, + 'patch_icpp(2)%alpha(1)': 0.1, + 'patch_icpp(3)%pres': 500000.0, + 'patch_icpp(3)%alpha_rho(1)': 900, + 'patch_icpp(3)%alpha(1)': 0.9, + 'fluid_pp(1)%gamma': 0.3, + 'fluid_pp(1)%pi_inf': 780000.0, +# 'Ca': 0.9769178386380458, +# 'Web': 13.927835051546392, +# 'Re_inv': 0.009954269975623245, + 'pref': 101325.0, + 'rhoref': 1000.0, +# 'bubble_model': 3, +# 'polytropic': 'T', +# 'polydisperse': 'F', +# 'thermal': 3, + 'R0ref': 1e-05, + 'patch_icpp(1)%r0': 1, + 'patch_icpp(1)%v0': 0, + 'patch_icpp(2)%r0': 1, + 'patch_icpp(2)%v0': 0, + 'patch_icpp(3)%r0': 1, + 'patch_icpp(3)%v0': 0, +# 'qbmm': 'F', +# 'dist_type': 2, +# 'poly_sigma': 0.3, +# 'R0_type': 1, +# 'sigR': 0.1, +# 'sigV': 0.1, +# 'rhoRV': 0.0, + 'x_domain%beg': 0.0, + 'x_domain%end': 1.0, + 'y_domain%beg': 0.0, + 'y_domain%end': 1.0, + 'bc_x%beg': -3, + 'bc_x%end': -3, + 'bc_y%beg': -3, + 'bc_y%end': -3, + 'patch_icpp(1)%geometry': 3, + 'patch_icpp(1)%y_centroid': 0.05, + 'patch_icpp(1)%length_y': 0.1, + 'patch_icpp(2)%y_centroid': 0.45, + 'patch_icpp(2)%length_y': 0.7, + 'patch_icpp(3)%y_centroid': 0.9, + 'patch_icpp(3)%length_y': 0.2, + 'patch_icpp(1)%x_centroid': 0.5, + 'patch_icpp(1)%length_x': 1, + 'patch_icpp(1)%vel(1)': 0.0, + 'patch_icpp(1)%vel(2)': 0.0, + 'patch_icpp(2)%geometry': 3, + 'patch_icpp(2)%x_centroid': 0.5, + 'patch_icpp(2)%length_x': 1, + 'patch_icpp(2)%vel(1)': 0.0, + 'patch_icpp(2)%vel(2)': 0.0, + 'patch_icpp(3)%geometry': 3, + 'patch_icpp(3)%x_centroid': 0.5, + 'patch_icpp(3)%length_x': 1, + 'patch_icpp(3)%vel(1)': 0.0, + 'patch_icpp(3)%vel(2)': 0.0, + 'hyperelasticity': 'T', + 'hyper_model': 1, + 'fd_order': 4, + 'patch_icpp(1)%tau_e(1)': 0.0, + 'patch_icpp(2)%tau_e(1)': 0.0, + 'patch_icpp(3)%tau_e(1)': 0.0, + 'fluid_pp(1)%G': 1000000000.0, + 'fluid_pp(2)%gamma': 0.3, + 'fluid_pp(2)%pi_inf': 780000.0, + 'patch_icpp(1)%alpha_rho(2)': 100, + 'patch_icpp(1)%alpha(2)': 0.1, + 'patch_icpp(2)%alpha_rho(2)': 900, + 'patch_icpp(2)%alpha(2)': 0.9, + 'patch_icpp(3)%alpha_rho(2)': 100, + 'patch_icpp(3)%alpha(2)': 0.1, + 'fluid_pp(2)%G': 0.0, #50000.0, + 'patch_icpp(1)%tau_e(2)': 0.0, + 'patch_icpp(1)%tau_e(3)': 0.0, + 'patch_icpp(2)%tau_e(2)': 0.0, + 'patch_icpp(2)%tau_e(3)': 0.0, + 'patch_icpp(3)%tau_e(2)': 0.0, + 'patch_icpp(3)%tau_e(3)': 0.0, + 'parallel_io' : 'T', 'cons_vars_wrt' : 'T', + 'prim_vars_wrt': 'T', 'alpha_rho_wrt(1)': 'T', + 'rho_wrt' : 'T', 'mom_wrt(1)' : 'T', + 'vel_wrt(1)' : 'T', 'E_wrt' : 'T', + 'pres_wrt' : 'T', 'alpha_wrt(1)' : 'T', + 'gamma_wrt' : 'T', 'heat_ratio_wrt' : 'T', + 'pi_inf_wrt' : 'T', 'pres_inf_wrt' : 'T', + 'c_wrt' : 'T', + })) + diff --git a/examples/3D_hyper_1material/case.py b/examples/3D_hyper_1material/case.py new file mode 100644 index 0000000000..ee365445bc --- /dev/null +++ b/examples/3D_hyper_1material/case.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python3 + +# 3D -> Hyperelasticity -> 1 Fluid(s) + +import math, json + +print(json.dumps({ + 'run_time_info': 'T', + 'm': 24, + 'n': 24, + 'p': 24, + 'dt': 1e-07, + 't_step_start': 0, + 't_step_stop': 50, + 't_step_save': 25, +# 'cfl_adap_dt': 'T', +# 'cfl_target': 0.25, +# 'n_start': 0, +# 't_stop': 1e-06, +# 't_save': 1e-07, + 'num_patches': 3, + 'model_eqns': 3, + 'alt_soundspeed': 'F', + 'num_fluids': 1, + 'mpp_lim': 'F', + 'mixture_err': 'F', + 'time_stepper': 3, + 'weno_order': 5, + 'weno_eps': 1e-16, + 'mapped_weno': 'F', + 'null_weights': 'F', + 'mp_weno': 'F', + 'riemann_solver': 2, + 'wave_speeds': 1, + 'avg_state': 2, + 'format': 1, + 'precision': 2, + 'prim_vars_wrt': 'F', + 'parallel_io': 'F', + 'patch_icpp(1)%pres': 1000000.0, + 'patch_icpp(1)%alpha_rho(1)': 1000.0, + 'patch_icpp(1)%alpha(1)': 1.0, + 'patch_icpp(2)%pres': 100000.0, + 'patch_icpp(2)%alpha_rho(1)': 1000.0, + 'patch_icpp(2)%alpha(1)': 1.0, + 'patch_icpp(3)%pres': 500000.0, + 'patch_icpp(3)%alpha_rho(1)': 1000.0, + 'patch_icpp(3)%alpha(1)': 1.0, + 'fluid_pp(1)%gamma': 0.3, + 'fluid_pp(1)%pi_inf': 780000.0, +# 'Ca': 0.9769178386380458, +# 'Web': 13.927835051546392, +# 'Re_inv': 0.009954269975623245, + 'pref': 101325.0, + 'rhoref': 1000.0, +# 'bubble_model': 3, +# 'polytropic': 'T', +# 'polydisperse': 'F', +# 'thermal': 3, + 'R0ref': 1e-05, + 'patch_icpp(1)%r0': 1, + 'patch_icpp(1)%v0': 0, + 'patch_icpp(2)%r0': 1, + 'patch_icpp(2)%v0': 0, + 'patch_icpp(3)%r0': 1, + 'patch_icpp(3)%v0': 0, +# 'qbmm': 'F', +# 'dist_type': 2, +# 'poly_sigma': 0.3, +# 'R0_type': 1, +# 'sigR': 0.1, +# 'sigV': 0.1, +# 'rhoRV': 0.0, + 'x_domain%beg': 0.0, + 'x_domain%end': 1.0, + 'y_domain%beg': 0.0, + 'y_domain%end': 1.0, + 'z_domain%beg': 0.0, + 'z_domain%end': 1.0, + 'bc_x%beg': -3, + 'bc_x%end': -3, + 'bc_y%beg': -3, + 'bc_y%end': -3, + 'bc_z%beg': -3, + 'bc_z%end': -3, + 'patch_icpp(1)%geometry': 9, + 'patch_icpp(1)%z_centroid': 0.05, + 'patch_icpp(1)%length_z': 0.1, + 'patch_icpp(2)%z_centroid': 0.45, + 'patch_icpp(2)%length_z': 0.7, + 'patch_icpp(3)%z_centroid': 0.9, + 'patch_icpp(3)%length_z': 0.2, + 'patch_icpp(1)%y_centroid': 0.5, + 'patch_icpp(1)%length_y': 1, + 'patch_icpp(1)%x_centroid': 0.5, + 'patch_icpp(1)%length_x': 1, + 'patch_icpp(1)%vel(1)': 0.0, + 'patch_icpp(1)%vel(2)': 0.0, + 'patch_icpp(1)%vel(3)': 0.0, + 'patch_icpp(2)%geometry': 9, + 'patch_icpp(2)%y_centroid': 0.5, + 'patch_icpp(2)%length_y': 1, + 'patch_icpp(2)%x_centroid': 0.5, + 'patch_icpp(2)%length_x': 1, + 'patch_icpp(2)%vel(1)': 0.0, + 'patch_icpp(2)%vel(2)': 0.0, + 'patch_icpp(2)%vel(3)': 0.0, + 'patch_icpp(3)%geometry': 9, + 'patch_icpp(3)%y_centroid': 0.5, + 'patch_icpp(3)%length_y': 1, + 'patch_icpp(3)%x_centroid': 0.5, + 'patch_icpp(3)%length_x': 1, + 'patch_icpp(3)%vel(1)': 0.0, + 'patch_icpp(3)%vel(2)': 0.0, + 'patch_icpp(3)%vel(3)': 0.0, + 'hyperelasticity': 'T', + 'hyper_model': 1, + 'fd_order': 4, + 'patch_icpp(1)%tau_e(1)': 0.0, + 'patch_icpp(2)%tau_e(1)': 0.0, + 'patch_icpp(3)%tau_e(1)': 0.0, + 'fluid_pp(1)%G': 1000000000.0, + 'patch_icpp(1)%tau_e(2)': 0.0, + 'patch_icpp(1)%tau_e(3)': 0.0, + 'patch_icpp(2)%tau_e(2)': 0.0, + 'patch_icpp(2)%tau_e(3)': 0.0, + 'patch_icpp(3)%tau_e(2)': 0.0, + 'patch_icpp(3)%tau_e(3)': 0.0, + 'patch_icpp(1)%tau_e(4)': 0.0, + 'patch_icpp(1)%tau_e(5)': 0.0, + 'patch_icpp(1)%tau_e(6)': 0.0, + 'patch_icpp(2)%tau_e(4)': 0.0, + 'patch_icpp(2)%tau_e(5)': 0.0, + 'patch_icpp(2)%tau_e(6)': 0.0, + 'patch_icpp(3)%tau_e(4)': 0.0, + 'patch_icpp(3)%tau_e(5)': 0.0, + 'patch_icpp(3)%tau_e(6)': 0.0, + 'parallel_io' : 'T', 'cons_vars_wrt' : 'T', + 'prim_vars_wrt': 'T', 'alpha_rho_wrt(1)': 'T', + 'rho_wrt' : 'T', 'mom_wrt(1)' : 'T', + 'vel_wrt(1)' : 'T', 'E_wrt' : 'T', + 'pres_wrt' : 'T', 'alpha_wrt(1)' : 'T', + 'gamma_wrt' : 'T', 'heat_ratio_wrt' : 'T', + 'pi_inf_wrt' : 'T', 'pres_inf_wrt' : 'T', + 'c_wrt' : 'T', +})) diff --git a/examples/3D_hyper_2materials/case.py b/examples/3D_hyper_2materials/case.py new file mode 100644 index 0000000000..33857671aa --- /dev/null +++ b/examples/3D_hyper_2materials/case.py @@ -0,0 +1,150 @@ +#!/usr/bin/env python3 + +# 3D -> Hyperelasticity -> 2 Fluid(s) + +import math, json + +print(json.dumps({ + 'run_time_info': 'T', + 'm': 24, + 'n': 24, + 'p': 24, + 'dt': 1e-06, + 't_step_start': 0, + 't_step_stop': 50, + 't_step_save': 25, + 'num_patches': 3, + 'model_eqns': 3, + 'alt_soundspeed': 'F', + 'num_fluids': 2, + 'mpp_lim': 'F', + 'mixture_err': 'F', + 'time_stepper': 3, + 'weno_order': 5, + 'weno_eps': 1e-16, + 'mapped_weno': 'F', + 'null_weights': 'F', + 'mp_weno': 'F', + 'riemann_solver': 2, + 'wave_speeds': 1, + 'avg_state': 2, + 'format': 1, + 'precision': 2, + 'prim_vars_wrt': 'F', + 'parallel_io': 'F', + 'patch_icpp(1)%pres': 1000000.0, + 'patch_icpp(1)%alpha_rho(1)': 900.0, + 'patch_icpp(1)%alpha(1)': 0.9, + 'patch_icpp(2)%pres': 100000.0, + 'patch_icpp(2)%alpha_rho(1)': 100, + 'patch_icpp(2)%alpha(1)': 0.1, + 'patch_icpp(3)%pres': 500000.0, + 'patch_icpp(3)%alpha_rho(1)': 900, + 'patch_icpp(3)%alpha(1)': 0.9, + 'fluid_pp(1)%gamma': 0.3, + 'fluid_pp(1)%pi_inf': 780000.0, +# 'Ca': 0.9769178386380458, +# 'Web': 13.927835051546392, +# 'Re_inv': 0.009954269975623245, + 'pref': 101325.0, + 'rhoref': 1000.0, +# 'bubble_model': 3, +# 'polytropic': 'T', +# 'polydisperse': 'F', +# 'thermal': 3, + 'R0ref': 1e-05, + 'patch_icpp(1)%r0': 1, + 'patch_icpp(1)%v0': 0, + 'patch_icpp(2)%r0': 1, + 'patch_icpp(2)%v0': 0, + 'patch_icpp(3)%r0': 1, + 'patch_icpp(3)%v0': 0, +# 'qbmm': 'F', +# 'dist_type': 2, +# 'poly_sigma': 0.3, +# 'R0_type': 1, +# 'sigR': 0.1, +# 'sigV': 0.1, +# 'rhoRV': 0.0, + 'x_domain%beg': 0.0, + 'x_domain%end': 1.0, + 'y_domain%beg': 0.0, + 'y_domain%end': 1.0, + 'z_domain%beg': 0.0, + 'z_domain%end': 1.0, + 'bc_x%beg': -3, + 'bc_x%end': -3, + 'bc_y%beg': -3, + 'bc_y%end': -3, + 'bc_z%beg': -3, + 'bc_z%end': -3, + 'patch_icpp(1)%geometry': 9, + 'patch_icpp(1)%z_centroid': 0.05, + 'patch_icpp(1)%length_z': 0.1, + 'patch_icpp(2)%z_centroid': 0.45, + 'patch_icpp(2)%length_z': 0.7, + 'patch_icpp(3)%z_centroid': 0.9, + 'patch_icpp(3)%length_z': 0.2, + 'patch_icpp(1)%y_centroid': 0.5, + 'patch_icpp(1)%length_y': 1, + 'patch_icpp(1)%x_centroid': 0.5, + 'patch_icpp(1)%length_x': 1, + 'patch_icpp(1)%vel(1)': 0.0, + 'patch_icpp(1)%vel(2)': 0.0, + 'patch_icpp(1)%vel(3)': 0.0, + 'patch_icpp(2)%geometry': 9, + 'patch_icpp(2)%y_centroid': 0.5, + 'patch_icpp(2)%length_y': 1, + 'patch_icpp(2)%x_centroid': 0.5, + 'patch_icpp(2)%length_x': 1, + 'patch_icpp(2)%vel(1)': 0.0, + 'patch_icpp(2)%vel(2)': 0.0, + 'patch_icpp(2)%vel(3)': 0.0, + 'patch_icpp(3)%geometry': 9, + 'patch_icpp(3)%y_centroid': 0.5, + 'patch_icpp(3)%length_y': 1, + 'patch_icpp(3)%x_centroid': 0.5, + 'patch_icpp(3)%length_x': 1, + 'patch_icpp(3)%vel(1)': 0.0, + 'patch_icpp(3)%vel(2)': 0.0, + 'patch_icpp(3)%vel(3)': 0.0, + 'hyperelasticity': 'T', + 'hyper_model': 1, + 'fd_order': 4, + 'patch_icpp(1)%tau_e(1)': 0.0, + 'patch_icpp(2)%tau_e(1)': 0.0, + 'patch_icpp(3)%tau_e(1)': 0.0, + 'fluid_pp(1)%G': 1000000000.0, + 'fluid_pp(2)%gamma': 0.3, + 'fluid_pp(2)%pi_inf': 780000.0, + 'patch_icpp(1)%alpha_rho(2)': 100, + 'patch_icpp(1)%alpha(2)': 0.1, + 'patch_icpp(2)%alpha_rho(2)': 900, + 'patch_icpp(2)%alpha(2)': 0.9, + 'patch_icpp(3)%alpha_rho(2)': 100, + 'patch_icpp(3)%alpha(2)': 0.1, + 'fluid_pp(2)%G': 0.0, #50000.0, + 'patch_icpp(1)%tau_e(2)': 0.0, + 'patch_icpp(1)%tau_e(3)': 0.0, + 'patch_icpp(2)%tau_e(2)': 0.0, + 'patch_icpp(2)%tau_e(3)': 0.0, + 'patch_icpp(3)%tau_e(2)': 0.0, + 'patch_icpp(3)%tau_e(3)': 0.0, + 'patch_icpp(1)%tau_e(4)': 0.0, + 'patch_icpp(1)%tau_e(5)': 0.0, + 'patch_icpp(1)%tau_e(6)': 0.0, + 'patch_icpp(2)%tau_e(4)': 0.0, + 'patch_icpp(2)%tau_e(5)': 0.0, + 'patch_icpp(2)%tau_e(6)': 0.0, + 'patch_icpp(3)%tau_e(4)': 0.0, + 'patch_icpp(3)%tau_e(5)': 0.0, + 'patch_icpp(3)%tau_e(6)': 0.0, + 'parallel_io' : 'T', 'cons_vars_wrt' : 'T', + 'prim_vars_wrt': 'T', 'alpha_rho_wrt(1)': 'T', + 'rho_wrt' : 'T', 'mom_wrt(1)' : 'T', + 'vel_wrt(1)' : 'T', 'E_wrt' : 'T', + 'pres_wrt' : 'T', 'alpha_wrt(1)' : 'T', + 'gamma_wrt' : 'T', 'heat_ratio_wrt' : 'T', + 'pi_inf_wrt' : 'T', 'pres_inf_wrt' : 'T', + 'c_wrt' : 'T', +})) 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_checker_common.fpp b/src/common/m_checker_common.fpp index 2b8c1a4d12..a0d4e41aa5 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 @@ -132,25 +132,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 @@ -177,6 +158,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 17d644259b..1379081baa 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -1039,6 +1039,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) @@ -1053,8 +1063,6 @@ contains end do !$acc end parallel loop - !print *, 'I got here AA' - end subroutine s_convert_conservative_to_primitive_variables !> The following procedure handles the conversion between 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 5faf15f597..33f226d0b1 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 @@ -343,7 +345,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 ! Volume fraction(s) if ((model_eqns == 2) .or. (model_eqns == 3)) then @@ -604,6 +606,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 @@ -1371,6 +1392,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 edb636ad2b..86855b5343 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -487,21 +487,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 17e1d450ea..a3b9a6ebd6 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -635,24 +635,6 @@ contains end if end if - if (hypoelasticity .or. hyperelasticity) then - elasticity = .true. - stress_idx%beg = sys_size + 1 - stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2 - ! number of 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 - 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 (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx @@ -675,24 +657,6 @@ contains internalEnergies_idx%end = adv_idx%end + num_fluids sys_size = internalEnergies_idx%end - if (hypoelasticity .or. hyperelasticity) then - elasticity = .true. - stress_idx%beg = sys_size + 1 - stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2 - ! number of 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 - 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 (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx @@ -756,6 +720,26 @@ contains end if end if + 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 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 + 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 @@ -852,11 +836,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_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index f5655bf025..5a16478dc0 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -55,7 +55,7 @@ contains & 'mixlayer_perturb', 'bubbles_euler', 'polytropic', 'polydisperse',& & 'qbmm', 'file_per_process', 'adv_n', 'ib' , 'cfl_adap_dt', & & 'cfl_const_dt', 'cfl_dt', 'surface_tension', & - & 'hyperelasticity', 'pre_stress' ] + & 'hyperelasticity', 'pre_stress'] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp index 5660761d63..b995bc59a2 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/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index a23519a37b..6e5c45bffb 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -140,7 +140,8 @@ contains file_per_process, relax, relax_model, & palpha_eps, ptgalpha_eps, ib, num_ibs, patch_ib, & sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, & - n_start_old, surface_tension, hyperelasticity, pre_stress, rkck_adap_dt + n_start_old, surface_tension, hyperelasticity, pre_stress, & + rkck_adap_dt ! Inquiring the status of the pre_process.inp file file_loc = 'pre_process.inp' diff --git a/src/simulation/m_boundary_conditions.fpp b/src/simulation/m_boundary_conditions.fpp index 073c62d072..b0b96d11b2 100644 --- a/src/simulation/m_boundary_conditions.fpp +++ b/src/simulation/m_boundary_conditions.fpp @@ -25,6 +25,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), dimension(startx:, starty:, startz:, 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 @@ -236,6 +238,26 @@ contains end do end do + 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 + else !< bc_x%end !$acc parallel loop collapse(4) gang vector default(present) @@ -250,6 +272,25 @@ contains end do end do + 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 end if !< y-direction @@ -269,6 +310,26 @@ contains end do end do + 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 + else !< bc_y%end !$acc parallel loop collapse(4) gang vector default(present) @@ -283,6 +344,26 @@ contains end do end do + 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 + end if !< z-direction @@ -302,6 +383,25 @@ contains end do end do + 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 else !< bc_z%end !$acc parallel loop collapse(4) gang vector default(present) @@ -316,6 +416,26 @@ contains end do end do + 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 + end if end if diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index 278e4c055e..6ec4ca746f 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 a7f6665df4..d49bb7fd98 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 @@ -122,7 +121,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") @@ -258,9 +257,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 c2e6ddba27..4a16d398f8 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 6963d6dfe2..3a86ee4efd 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -690,7 +690,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. @@ -891,26 +891,6 @@ contains end if end if - if (hypoelasticity .or. hyperelasticity) then - elasticity = .true. - 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 - hyper_model = 1 - end if - if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx @@ -929,25 +909,6 @@ contains internalEnergies_idx%end = adv_idx%end + num_fluids sys_size = internalEnergies_idx%end - if (hypoelasticity .or. hyperelasticity) then - elasticity = .true. - stress_idx%beg = sys_size + 1 - stress_idx%end = sys_size + (num_dims*(num_dims + 1))/2 - ! number of 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 (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx @@ -1045,6 +1006,27 @@ contains end if ! 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 @@ -1235,11 +1217,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..1c6bca060f 100644 --- a/src/simulation/m_hyperelastic.fpp +++ b/src/simulation/m_hyperelastic.fpp @@ -116,103 +116,158 @@ 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 + + !$acc loop seq + do i = 1, tensor_size + tensora(i) = 0._wp + end do if (G > verysmall) then - !$acc loop seq - do i = 1, tensor_size - tensora(i) = 0_wp - 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 + ! STEP 1: computing grad_xi (tensora) using finite differences + ! tensora(1,2): dxix_x, dxiy_x + ! tensora(3,4): dxix_dy, dxiy_dy + ! 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 (p > 0) then + ! 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) + elseif (n > 0) then + ! 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) + ! 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) + ! print *, 'tensora(1)::', tensora,'tensora(2) ::', tensora(2), 'tensora(3) ::', tensora(3), 'tensora(4)::', tensora(4) + else + ! derivatives in the x-direction + tensora(1) = tensora(1) + q_prim_vf(xibeg)%sf(j + r, k, l)*fd_coeff_x(r, j) + 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)) - - if (tensorb(tensor_size) > verysmall) then - ! STEP 2c: computing the inverse of grad_xi tensor = F - ! tensorb is the adjoint, tensora becomes F - !$acc loop seq - do i = 1, tensor_size - 1 - tensora(i) = tensorb(i)/tensorb(tensor_size) - 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 - btensor%vf(b_size)%sf(j, k, l) = tensorb(tensor_size) - ! STEP 5a: 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) + + 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 + tensora(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)) + if (tensora(tensor_size) > verysmall) then + ! STEP 2c: computing the inverse of grad_xi tensor = F (tensora) + !$acc loop seq + do i = 1, tensor_size - 1 + tensora(i) = tensorb(i)/tensora(tensor_size) + end do + ! STEP 2d: computing J = det(F) + tensorb(tensor_size) = 1._wp/tensora(tensor_size) + ! 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 + elseif (n > 0) then + ! STEP 2a: computing the cofactor (tensorb) of the grad_xi tensor for the inverse + tensorb(1) = tensora(4) + tensorb(2) = -tensora(2) + tensorb(3) = -tensora(3) + tensorb(4) = tensora(1) + ! STEP 2b: computing the determinant of the grad_xi tensor + tensora(tensor_size) = tensora(1)*tensora(4) - tensora(3)*tensora(2) + ! print *, 'I compute determinant::', tensora(tensor_size) + if (tensora(tensor_size) > verysmall) then + ! STEP 2c: computing the inverse of grad_xi tensor = F (tensora) + !$acc loop seq + do i = 1, tensor_size - 1 + tensora(i) = tensorb(i)/tensora(tensor_size) + end do + ! STEP 2d: computing J = det(F) + tensorb(tensor_size) = 1._wp/tensora(tensor_size) + ! print *, 'I compute J::', tensorb(tensor_size) + ! STEP 2e: override adjoint (tensorb) to be F transpose F + tensorb(1) = tensora(4)**2 + tensora(3)**2 + tensorb(4) = tensora(2)**2 + tensora(1)**2 + tensorb(2) = (-tensora(2))*tensora(4) + tensora(1)*(-tensora(3)) + tensorb(3) = tensorb(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 - ! STEP 5b: 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 - !$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) - end do + else + ! 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) + if (tensora(tensor_size) > verysmall) then + ! STEP 3: update the btensor b_xx, this is consistent with Riemann solvers + btensor%vf(1)%sf(j, k, l) = tensorb(1)**2 + end if + 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 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 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 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) + end do end if end do 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,21 +283,32 @@ contains real(wp) :: trace real(wp), parameter :: f13 = 1._wp/3._wp - integer :: i !< Generic loop iterators - - ! 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) + integer :: i + + 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 + elseif (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 + else + ! calculate the deviatoric of the tensor + btensor(1)%sf(j, k, l) = 2._wp*btensor(1)%sf(j, k, l)/3._wp + end if - ! 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 ! dividing by the jacobian for neo-Hookean model ! setting the tensor to the stresses for riemann solver !$acc loop seq do i = 1, b_size - 1 q_prim_vf(strxb + i - 1)%sf(j, k, l) = & - G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) + G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) end do ! compute the invariant without the elastic modulus q_prim_vf(xiend + 1)%sf(j, k, l) = & @@ -250,7 +316,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,23 +334,32 @@ contains real(wp) :: trace real(wp), parameter :: f13 = 1._wp/3._wp - integer :: i !< Generic loop iterators - - !TODO Make this 1D and 2D capable - ! 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 - 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 - btensor(6)%sf(j, k, l) = btensor(6)%sf(j, k, l) - f13*trace + integer :: i + + 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 + elseif (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 + else + ! calculate the deviatoric of the tensor + btensor(1)%sf(j, k, l) = 2._wp*btensor(1)%sf(j, k, l)/3._wp + end if ! dividing by the jacobian for neo-Hookean model ! setting the tensor to the stresses for riemann solver !$acc loop seq do i = 1, b_size - 1 q_prim_vf(strxb + i - 1)%sf(j, k, l) = & - G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) + G*btensor(i)%sf(j, k, l)/btensor(b_size)%sf(j, k, l) end do ! compute the invariant without the elastic modulus q_prim_vf(xiend + 1)%sf(j, k, l) = & diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp index 2e676ba087..e7a7bbd059 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 1764c06832..f0450f13b6 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -179,7 +179,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 d9617f9304..b00f4a0b06 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -666,7 +666,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 @@ -705,7 +708,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 bbb74de381..65d3d06e28 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -465,7 +465,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 @@ -544,12 +543,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) ! Additional terms in 2D and 3D @@ -561,37 +561,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 ! Enthalpy with elastic energy H_L = (E_L + pres_L)/rho_L @@ -648,6 +642,7 @@ 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)) - & @@ -705,7 +700,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_dims flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(i)) = & @@ -744,7 +739,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) = & @@ -785,9 +780,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)) & @@ -799,7 +794,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) = & @@ -812,18 +819,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_dims @@ -1168,7 +1163,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 @@ -1190,10 +1184,9 @@ contains 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 - ! ENERGY ADJUSTMENTS FOR HYPOELASTIC ENERGY + ! energy adjustments for hypoelastic energy if (hypoelasticity) then !$acc loop seq do i = 1, strxe - strxb + 1 @@ -1221,7 +1214,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 @@ -1270,7 +1263,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 + & @@ -1345,8 +1338,9 @@ contains 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) = & @@ -1354,7 +1348,7 @@ 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 @@ -1363,21 +1357,21 @@ contains (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; !$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)))))) - & @@ -1387,7 +1381,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) = & @@ -1395,7 +1389,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) @@ -1404,7 +1398,7 @@ 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 @@ -1422,7 +1416,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 @@ -1432,7 +1426,7 @@ contains end do end if - ! REFERENCE MAP FLUX. + ! reference map flux. if (hyperelasticity) then !$acc loop seq do i = 1, num_dims @@ -1444,7 +1438,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) + & @@ -1485,7 +1479,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 @@ -1845,7 +1838,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 @@ -1853,14 +1845,12 @@ 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 H_L = (E_L + pres_L)/rho_L @@ -2382,7 +2372,7 @@ contains 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 @@ -2411,7 +2401,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 @@ -2521,8 +2511,9 @@ contains xi_M = (5e-1_wp + sign(5e-1_wp, s_S)) xi_P = (5e-1_wp - sign(5e-1_wp, s_S)) - ! COMPUTING THE HLLC FLUXES - ! MASS FLUX. + ! Computing the hllc fluxes + + ! mass flux. if (low_Mach == 1) then @:compute_low_Mach_correction() else @@ -2538,7 +2529,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 @@ -2559,7 +2550,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) + & @@ -2572,17 +2563,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)))))) - & @@ -2592,7 +2585,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 @@ -2602,7 +2595,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) = & @@ -2612,7 +2605,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) @@ -2625,7 +2618,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 diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index 0c6f95a672..2fc787cf5b 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -77,7 +77,7 @@ subroutine s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, E = gamma*pres + pi_inf + 5e-1_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 diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 0a84b34118..f3c78215a7 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -167,11 +167,11 @@ 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 + rkck_adap_dt, rkck_tolerance ! Checking that an input file has been provided by the user. If it ! has, then the input file is read in, otherwise, simulation exits. diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 2ffbc8e58b..40d907e955 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -843,12 +843,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) @@ -1067,7 +1067,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) @@ -1081,7 +1081,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) @@ -1095,7 +1095,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) @@ -1109,7 +1109,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) @@ -1123,7 +1123,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) @@ -1137,7 +1137,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) @@ -1154,7 +1154,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 d037f39af0..c695e08385 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -244,6 +244,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, }) @@ -362,6 +363,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, @@ -383,6 +385,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 8eb3819a84..8015f989e2 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -526,7 +526,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, @@ -848,8 +848,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