Skip to content

Commit

Permalink
Add Body Forces (#449)
Browse files Browse the repository at this point in the history 
Co-authored-by: Ben Wilfong <bwilfong3@gatech.edu>
  • Loading branch information
@wilfonba
wilfonba and wilfonba authored Jun 3, 2024
1 parent 9369312 commit daa8e85
Show file tree
Hide file tree
Showing 26 changed files with 1,235 additions and 12 deletions.
18 changes: 17 additions & 1 deletion docs/documentation/case.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -668,8 +668,24 @@ The parameters are optionally used to define initial velocity profiles and pertu
| ---: | :----: | :--- |
| `pi_fac` | Real | Ratio of artificial and true `pi_\infty` values|

- `pi_fac` specifies the ratio of artificial and true `pi_\infty` values (`=` artificial `pi_\infty` / true `pi_\infty`). This parameter enables the use of true `pi_\infty` in bubble dynamics models, when the `pi_\infty` given in the `case.py` file is an artificial value.
- `pi_fac` specifies the ratio of artificial and true `pi_\infty` values (`=` artificial `pi_\infty` / true `pi_\infty`). i
This parameter enables the use of true `pi_\infty` in bubble dynamics models, when the `pi_\infty` given in the `case.py` file is an artificial value.

### 13. Body Forces

| Parameter | Type | Description |
| ---: | :---: | :--- |
| `bf_x[y,z]` | Logical | Enable body forces in the x[y,z] direction |
| `k_x[y,y]` | Real | Magnitude of oscillating acceleration |
| `w_x[y,z]` | Real | Frequency of oscillating acceleration |
| `p_x[y,z]` | Real | Phase shift of oscillating acceleration |
| `g_x[y,z]` | Real | Magnitude of bacground acceleration |

`k_x[y,z]`, `w_x[y,z]`, `p_x[y,z]`, and `g_x[y,z]` define an oscillating acceleration in the `x[y,z]` direction with the form

$$ a_{x[y,z]} = g_{x[y,z]} + k_{x[y,z]}\sin\left(w_{x[y,z]}t + p_{x[y,z]}\right). $$

Positive accelerations are in the `x[y,z]` direction are in the positive `x[y,z]` direction by convention.

## Enumerations

Expand Down
9 changes: 9 additions & 0 deletions examples/2D_rayleigh_taylor/README.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
# Rayleigh-Taylor Instability (2D)

## Final Condition

![Final Condition](final_condition.png)

## Centerline Velocities

![Linear Theory Comparision](linear_theory.jpg)

Check warning on line 9 in examples/2D_rayleigh_taylor/README.md

View workflow job for this annotation

GitHub Actions / Spell Check 

"Comparision" should be "Comparison".
121 changes: 121 additions & 0 deletions examples/2D_rayleigh_taylor/case.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,121 @@
#!/usr/bin/env python3

import math
import json

lam = 0.2
h = 1.2
k = 2*math.pi/lam
amp = 0.05/k

# Numerical setup
x0 = -lam/2
x1 = lam/2
y0 = 0.
y1 = h

Nx = 199
Ny = 1199

eps = 1e-6

dx = lam/(Nx + 1)
c = math.sqrt(1.4*1e5/1)
cfl = 0.5
dt = cfl*dx/c

Nt = math.ceil(2/dt)
Ns = math.ceil(Nt/100)

#Configuration case dictionary
data = {
# Logistics =============================
#'case_dir' : '\'.\'',
'run_time_info' : 'T',
# =======================================

# Computational Domain ==================
'x_domain%beg' : x0,
'x_domain%end' : x1,
'y_domain%beg' : y0,
'y_domain%end' : y1,
'm' : Nx,
'n' : Ny,
'p' : 0,
'cyl_coord' : 'F',
'dt' : dt,
't_step_start' : 0,
't_step_stop' : Nt,
't_step_save' : Ns,
# =======================================

# Simulation Algorithm ==================
'model_eqns' : 3,
'alt_soundspeed' : 'F',
'adv_alphan' : 'T',
'mixture_err' : 'T',
'mpp_lim' : 'T',
'time_stepper' : 3,
'avg_state' : 2,
'weno_order' : 5,
'weno_eps' : 1e-16,
'mapped_weno' : 'T',
'null_weights' : 'F',
'mp_weno' : 'T',
'weno_Re_flux' : 'T',
'riemann_solver' : 2,
'wave_speeds' : 1,
'bc_x%beg' : -2,
'bc_x%end' : -2,
'bc_y%beg' : -16,
'bc_y%end' : -16,
'num_patches' : 1,
'num_fluids' : 2,
# =======================================

# Database Structure Parameters =========
'format' : 1,
'precision' : 2,
'prim_vars_wrt' : 'T',
'parallel_io' : 'T',
# =======================================

# Fluid Parameters (Heavy Gas) ==========
'fluid_pp(1)%gamma' : 1.E+00/(1.4E+00-1.E+00),
'fluid_pp(1)%pi_inf' : 0.E+00,
'fluid_pp(1)%Re(1)' : 1/0.0219,
# =======================================

# Fluid Parameters (Light Gas) ==========
'fluid_pp(2)%gamma' : 1.E+00/(1.4E+00-1.E+00),
'fluid_pp(2)%pi_inf' : 0.E+00,
'fluid_pp(2)%Re(1)' : 1/0.0073,
# =======================================

# Body Forces ===========================
'bf_y' : 'T',
'k_y' : 0.,
'w_y' : 0.,
'p_y' : 0.,
'g_y' : -9.81,
# ======================================

# Water Patch ==========================
'patch_icpp(1)%geometry' : 7,
'patch_icpp(1)%hcid' : 204,
'patch_icpp(1)%x_centroid' : 0,
'patch_icpp(1)%y_centroid' : h/2,
'patch_icpp(1)%length_x' : lam,
'patch_icpp(1)%length_y' : h,
'patch_icpp(1)%vel(1)' : 0.0,
'patch_icpp(1)%vel(2)' : 0.0,
'patch_icpp(1)%pres' : 1e5,
'patch_icpp(1)%alpha_rho(1)': (1-eps),
'patch_icpp(1)%alpha_rho(2)': eps*1,
'patch_icpp(1)%alpha(1)' : 1-eps,
'patch_icpp(1)%alpha(2)' : eps,
# ======================================

}

print(json.dumps(data))
Binary file added examples/2D_rayleigh_taylor/final_condition.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added examples/2D_rayleigh_taylor/linear_theory.jpg
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
9 changes: 9 additions & 0 deletions examples/3D_rayleigh_taylor/README.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
# Rayleigh-Taylor Instability (3D)

## Final Condition

![Final Condition](final_condition.png)

## Centerline Velocities

![Linear Theory Comparision](linear_theory.jpg)

Check warning on line 9 in examples/3D_rayleigh_taylor/README.md

View workflow job for this annotation

GitHub Actions / Spell Check 

"Comparision" should be "Comparison".
131 changes: 131 additions & 0 deletions examples/3D_rayleigh_taylor/case.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,131 @@
#!/usr/bin/env python3

import math
import json

lam = 0.2
h = 1.2
k = 2*math.pi/lam
amp = 0.05/k

# Numerical setup
x0 = 0
x1 = lam/2
y0 = 0.
y1 = h
z0 = 0
z1 = lam/2

Nx = 99
Ny = 1199
Nz = 99

eps = 1e-6

dx = lam/2/(Nx + 1)
c = math.sqrt(1.4*1e5/1)
cfl = 0.4
dt = cfl*dx/c

Nt = math.ceil(0.2/dt)
Ns = math.ceil(Nt/100)

#Configuration case dictionary
data = {
# Logistics =============================
#'case_dir' : '\'.\'',
'run_time_info' : 'T',
# =======================================

# Computational Domain ==================
'x_domain%beg' : x0,
'x_domain%end' : x1,
'y_domain%beg' : y0,
'y_domain%end' : y1,
'z_domain%beg' : z0,
'z_domain%end' : z1,
'm' : Nx,
'n' : Ny,
'p' : Nz,
'cyl_coord' : 'F',
'dt' : dt,
't_step_start' : 0,
't_step_stop' : Nt,
't_step_save' : Ns,
# =======================================

# Simulation Algorithm ==================
'model_eqns' : 2,
'alt_soundspeed' : 'F',
'adv_alphan' : 'T',
'mixture_err' : 'T',
'mpp_lim' : 'T',
'time_stepper' : 3,
'avg_state' : 2,
'weno_order' : 5,
'weno_eps' : 1e-16,
'mapped_weno' : 'T',
'null_weights' : 'F',
'mp_weno' : 'T',
'weno_Re_flux' : 'T',
'riemann_solver' : 2,
'wave_speeds' : 1,
'bc_x%beg' : -2,
'bc_x%end' : -3,
'bc_y%beg' : -16,
'bc_y%end' : -16,
'bc_z%beg' : -2,
'bc_z%end' : -3,
'num_patches' : 1,
'num_fluids' : 2,
# =======================================

# Database Structure Parameters =========
'format' : 1,
'precision' : 2,
'prim_vars_wrt' : 'T',
'parallel_io' : 'T',
# =======================================

# Fluid Parameters (Heavy Gas) ==========
'fluid_pp(1)%gamma' : 1.E+00/(1.4E+00-1.E+00),
'fluid_pp(1)%pi_inf' : 0.E+00,
'fluid_pp(1)%Re(1)' : 1/0.0219,
# =======================================

# Fluid Parameters (Light Gas) ==========
'fluid_pp(2)%gamma' : 1.E+00/(1.4E+00-1.E+00),
'fluid_pp(2)%pi_inf' : 0.E+00,
'fluid_pp(2)%Re(1)' : 1/0.0073,
# =======================================

# Body Forces ===========================
'bf_y' : 'T',
'k_y' : 0.,
'w_y' : 0.,
'p_y' : 0.,
'g_y' : -98.1,
# ======================================

# Water Patch ==========================
'patch_icpp(1)%geometry' : 13,
'patch_icpp(1)%hcid' : 300,
'patch_icpp(1)%x_centroid' : 0,
'patch_icpp(1)%y_centroid' : h/2,
'patch_icpp(1)%z_centroid' : 0,
'patch_icpp(1)%length_x' : lam,
'patch_icpp(1)%length_y' : h,
'patch_icpp(1)%length_z' : h,
'patch_icpp(1)%vel(1)' : 0.0,
'patch_icpp(1)%vel(2)' : 0.0,
'patch_icpp(1)%vel(3)' : 0.0,
'patch_icpp(1)%pres' : 1e5,
'patch_icpp(1)%alpha_rho(1)': (1-eps),
'patch_icpp(1)%alpha_rho(2)': eps*1,
'patch_icpp(1)%alpha(1)' : 1-eps,
'patch_icpp(1)%alpha(2)' : eps,
# ======================================

}

print(json.dumps(data))
Binary file added examples/3D_rayleigh_taylor/final_condition.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added examples/3D_rayleigh_taylor/linear_theory.jpg
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
34 changes: 34 additions & 0 deletions src/pre_process/include/2dHardcodedIC.fpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@
real(kind(0d0)) :: eps
real(kind(0d0)) :: r, rmax, gam, umax, p0

real(kind(0d0)) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph

eps = 1e-9

#:enddef
Expand Down Expand Up @@ -66,6 +68,38 @@

q_prim_vf(contxb)%sf(i, j, 0) = q_prim_vf(E_idx)%sf(i, j, 0)**(1d0/gam)

case (204) ! Rayleigh-Taylor instability
rhoH = 3
rhoL = 1
pRef = 1e5
pInt = pRef
h = 0.7
lam = 0.2
wl = 2*pi/lam
amp = 0.05/wl

intH = amp*sin(2*pi*x_cc(i)/lam - pi/2) + h

alph = 5d-1*(1 + tanh((y_cc(j) - intH)/2.5e-3))

if (alph < eps) alph = eps
if (alph > 1 - eps) alph = 1 - eps

if (y_cc(j) > intH) then
q_prim_vf(advxb)%sf(i, j, 0) = alph
q_prim_vf(advxe)%sf(i, j, 0) = 1 - alph
q_prim_vf(contxb)%sf(i, j, 0) = alph*rhoH
q_prim_vf(contxe)%sf(i, j, 0) = (1 - alph)*rhoL
q_prim_vf(E_idx)%sf(i, j, 0) = pref + rhoH*9.81*(1.2 - y_cc(j))
else
q_prim_vf(advxb)%sf(i, j, 0) = alph
q_prim_vf(advxe)%sf(i, j, 0) = 1 - alph
q_prim_vf(contxb)%sf(i, j, 0) = alph*rhoH
q_prim_vf(contxe)%sf(i, j, 0) = (1 - alph)*rhoL
pInt = pref + rhoH*9.81*(1.2 - intH)
q_prim_vf(E_idx)%sf(i, j, 0) = pInt + rhoL*9.81*(intH - y_cc(j))
end if

case default
if (proc_rank == 0) then
call s_int_to_str(patch_id, iStr)
Expand Down
38 changes: 37 additions & 1 deletion src/pre_process/include/3dHardcodedIC.fpp
View file
Original file line number Diff line number Diff line change
@@ -1,12 +1,48 @@
#:def Hardcoded3DVariables()
! Place any declaration of intermediate variables here

real(kind(0d0)) :: rhoH, rhoL, pRef, pInt, h, lam, wl, amp, intH, alph

Check warning on line 4 in src/pre_process/include/3dHardcodedIC.fpp

View workflow job for this annotation

GitHub Actions / Spell Check 

"alph" should be "alpha".

real(kind(0d0)) :: eps

eps = 1e-9
#:enddef

#:def Hardcoded3D()

select case (patch_icpp(patch_id)%hcid)
case (300)
case (300) ! Rayleigh-Taylor instability
rhoH = 3
rhoL = 1
pRef = 1e5
pInt = pRef
h = 0.7
lam = 0.2
wl = 2*pi/lam
amp = 0.025/wl

intH = amp*(sin(2*pi*x_cc(i)/lam - pi/2) + sin(2*pi*z_cc(k)/lam - pi/2)) + h

alph = 5d-1*(1 + tanh((y_cc(j) - intH)/2.5e-3))

Check warning on line 26 in src/pre_process/include/3dHardcodedIC.fpp

View workflow job for this annotation

GitHub Actions / Spell Check 

"alph" should be "alpha".

if (alph < eps) alph = eps

Check warning on line 28 in src/pre_process/include/3dHardcodedIC.fpp

View workflow job for this annotation

GitHub Actions / Spell Check 

"alph" should be "alpha".

Check warning on line 28 in src/pre_process/include/3dHardcodedIC.fpp

View workflow job for this annotation

GitHub Actions / Spell Check 

"alph" should be "alpha".
if (alph > 1 - eps) alph = 1 - eps

Check warning on line 29 in src/pre_process/include/3dHardcodedIC.fpp

View workflow job for this annotation

GitHub Actions / Spell Check 

"alph" should be "alpha".

Check warning on line 29 in src/pre_process/include/3dHardcodedIC.fpp

View workflow job for this annotation

GitHub Actions / Spell Check 

"alph" should be "alpha".

if (y_cc(j) > intH) then
q_prim_vf(advxb)%sf(i, j, k) = alph

Check warning on line 32 in src/pre_process/include/3dHardcodedIC.fpp

View workflow job for this annotation

GitHub Actions / Spell Check 

"alph" should be "alpha".
q_prim_vf(advxe)%sf(i, j, k) = 1 - alph
q_prim_vf(contxb)%sf(i, j, k) = alph*rhoH
q_prim_vf(contxe)%sf(i, j, k) = (1 - alph)*rhoL
q_prim_vf(E_idx)%sf(i, j, k) = pref + rhoH*9.81*(1.2 - y_cc(j))
else
q_prim_vf(advxb)%sf(i, j, k) = alph
q_prim_vf(advxe)%sf(i, j, k) = 1 - alph
q_prim_vf(contxb)%sf(i, j, k) = alph*rhoH
q_prim_vf(contxe)%sf(i, j, k) = (1 - alph)*rhoL
pInt = pref + rhoH*9.81*(1.2 - intH)
q_prim_vf(E_idx)%sf(i, j, k) = pInt + rhoL*9.81*(intH - y_cc(j))
end if

! Put your variable assignments here
case default
call s_int_to_str(patch_id, iStr)
Expand Down
Loading

0 comments on commit daa8e85

Please sign in to comment.