Skip to content

Commit

Permalink
Merge pull request #205 from RasmitDevkota/master
Browse files Browse the repository at this point in the history
  • Loading branch information
sbryngelson authored Oct 2, 2023
2 parents 38ef7b9 + 0808431 commit 61f84b3
Show file tree
Hide file tree
Showing 10 changed files with 180 additions and 214 deletions.
2 changes: 1 addition & 1 deletion docs/documentation/case.md
Original file line number Diff line number Diff line change
Expand Up @@ -560,7 +560,7 @@ The entries labeled "Characteristic." are characteristic boundary conditions bas
| 3 | Rectangle | 2 | N | Coordinate-aligned. Requires `[x,y]_centroid` and `[x,y]_length`. |
| 4 | Sweep line | 2 | Y | Not coordinate aligned. Requires `[x,y]_centroid` and `normal(i)`. |
| 5 | Ellipse | 2 | Y | Requires `[x,y]_centroid` and `radii(i)`. |
| 6 | Vortex | 2 | N | Isentropic flow disturbance. Requires `[x,y]_centroid` and `radius`. |
| 6 | N/A | 2 | N | No longer exists. Empty. |
| 7 | 2D analytical | 2 | N | Assigns the primitive variables as analytical functions. |
| 8 | Sphere | 3 | Y | Requires `[x,y,z]_centroid` and `radius` |
| 9 | Cuboid | 3 | N | Coordinate-aligned. Requires `[x,y,z]_centroid` and `[x,y,z]_length`. |
Expand Down
11 changes: 11 additions & 0 deletions examples/2D_isentropicvortex/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
# Isentropic vortex problem (2D)

Reference: Coralic, V., & Colonius, T. (2014). Finite-volume Weno scheme for viscous compressible multicomponent flows. Journal of Computational Physics, 274, 95–121. https://doi.org/10.1016/j.jcp.2014.06.003

## Density

![Density](alpha_rho1.png)

## Density Norms

![Density Norms](density_norms.png)
Binary file added examples/2D_isentropicvortex/alpha_rho1.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
122 changes: 122 additions & 0 deletions examples/2D_isentropicvortex/case.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,122 @@
import math
import json

# Parameters
epsilon = '5d0'
alpha = '1d0'
gamma = '1.4'

# Initial conditions
vel1_i = '0d0'
vel2_i = '0d0'
T_i = '1d0'
alpha_rho1_i = '1d0'
pres_i = '1d0'

# Perturbations
vel1 = f'{vel1_i} + (y - yc)*({epsilon}/(2d0*pi))*' + \
f'exp({alpha}*(1d0 - (x - xc)**2d0 - (y - yc)**2d0))'
vel2 = f'{vel2_i} - (x - xc)*({epsilon}/(2d0*pi))*' + \
f'exp({alpha}*(1d0 - (x - xc)**2d0 - (y - yc)**2d0))'
alpha_rho1 = f'{alpha_rho1_i}*(1d0 - ({alpha_rho1_i}/{pres_i})*({epsilon}/(2d0*pi))*' + \
f'({epsilon}/(8d0*{alpha}*({gamma} + 1d0)*pi))*' + \
f'exp(2d0*{alpha}*(1d0 - (x - xc)**2d0' + \
f'- (y - yc)**2d0))' + \
f')**{gamma}'
pres = f'{pres_i}*(1d0 - ({alpha_rho1_i}/{pres_i})*({epsilon}/(2d0*pi))*' + \
f'({epsilon}/(8d0*{alpha}*({gamma} + 1d0)*pi))*' + \
f'exp(2d0*{alpha}*(1d0 - (x - xc)**2d0' + \
f'- (y - yc)**2d0))' + \
f')**({gamma} + 1d0)'

# Numerical setup
Nx = 399
dx = 10./(1.*(Nx+1))

c = 1.4**0.5
C = 0.3
mydt = C * dx / c * 0.01
Nt = 1/mydt

# Configuring case dictionary
print(json.dumps({
# Logistics ================================================================
'run_time_info' : 'T',
# ==========================================================================

# Computational Domain Parameters ==========================================
'x_domain%beg' : -3,
'x_domain%end' : 3,
'y_domain%beg' : -3,
'y_domain%end' : 3,
'stretch_x' : True,
'stretch_y' : True,
'loops_x' : 2,
'loops_y' : 2,
'a_x' : 1.03,
'a_y' : 1.03,
'x_a' : -1.5,
'y_a' : -1.5,
'x_b' : 1.5,
'y_b' : 1.5,
'm' : Nx,
'n' : Nx,
'p' : 0,
'dt' : mydt,
't_step_start' : 0,
't_step_stop' : 10,
't_step_save' : 1,
# ==========================================================================

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

# Formatted Database Files Structure Parameters ============================
'format' : 1,
'precision' : 2,
'prim_vars_wrt' :'T',
'parallel_io' :'F',
'omega_wrt(3)' :'T',
'fd_order' : 2,
# ==========================================================================

# Patch 1 ==================================================================
'patch_icpp(1)%geometry' : 3,
'patch_icpp(1)%x_centroid' : 0,
'patch_icpp(1)%y_centroid' : 0,
'patch_icpp(1)%length_x' : 10.,
'patch_icpp(1)%length_y' : 10.,
'patch_icpp(1)%vel(1)' : vel1,
'patch_icpp(1)%vel(2)' : vel2,
'patch_icpp(1)%pres' : pres,
'patch_icpp(1)%alpha_rho(1)' : alpha_rho1,
'patch_icpp(1)%alpha(1)' : 1.,
# ==========================================================================

# Fluids Physical Parameters ===============================================
'fluid_pp(1)%gamma' : 1.E+00/(1.4-1.E+00),
'fluid_pp(1)%pi_inf' : 0.0,
# ==========================================================================
}))

# ==============================================================================
Binary file added examples/2D_isentropicvortex/density_norms.png
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
3 changes: 1 addition & 2 deletions src/common/m_derived_types.fpp
Original file line number Diff line number Diff line change
Expand Up @@ -123,8 +123,7 @@ module m_derived_types
type(ic_model_parameters) :: model !< Model parameters

real(kind(0d0)) :: epsilon, beta !<
!! The isentropic vortex parameters administrating, respectively, both
!! the amplitude of the disturbance as well as its domain of influence.
!! The spherical harmonics eccentricity parameters.

real(kind(0d0)), dimension(3) :: normal !<
!! Normal vector indicating the orientation of the patch. It is specified
Expand Down
155 changes: 31 additions & 124 deletions src/pre_process/m_assign_variables.f90
Original file line number Diff line number Diff line change
Expand Up @@ -112,94 +112,42 @@ subroutine s_assign_patch_mixture_primitive_variables(patch_id, j, k, l, &
integer :: i !< generic loop operator

! Assigning the mixture primitive variables of a uniform state patch
if (patch_icpp(patch_id)%geometry /= 6) then

! Transferring the identity of the smoothing patch
smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id
! Transferring the identity of the smoothing patch
smooth_patch_id = patch_icpp(patch_id)%smooth_patch_id

! Density
q_prim_vf(1)%sf(j, k, l) = &
eta*patch_icpp(patch_id)%rho &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%rho

! Velocity
do i = 1, E_idx - mom_idx%beg
q_prim_vf(i + 1)%sf(j, k, l) = &
1d0/q_prim_vf(1)%sf(j, k, l)* &
(eta*patch_icpp(patch_id)%rho &
*patch_icpp(patch_id)%vel(i) &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%rho &
*patch_icpp(smooth_patch_id)%vel(i))
end do
! Density
q_prim_vf(1)%sf(j, k, l) = &
eta*patch_icpp(patch_id)%rho &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%rho

! Specific heat ratio function
q_prim_vf(gamma_idx)%sf(j, k, l) = &
eta*patch_icpp(patch_id)%gamma &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%gamma

! Pressure
q_prim_vf(E_idx)%sf(j, k, l) = &
1d0/q_prim_vf(gamma_idx)%sf(j, k, l)* &
(eta*patch_icpp(patch_id)%gamma &
*patch_icpp(patch_id)%pres &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%gamma &
*patch_icpp(smooth_patch_id)%pres)

! Liquid stiffness function
q_prim_vf(pi_inf_idx)%sf(j, k, l) = &
eta*patch_icpp(patch_id)%pi_inf &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%pi_inf

! Assigning mixture primitive variables of isentropic vortex patch
else
! Velocity
do i = 1, E_idx - mom_idx%beg
q_prim_vf(i + 1)%sf(j, k, l) = &
1d0/q_prim_vf(1)%sf(j, k, l)* &
(eta*patch_icpp(patch_id)%rho &
*patch_icpp(patch_id)%vel(i) &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%rho &
*patch_icpp(smooth_patch_id)%vel(i))
end do

! Transferring the centroid of the isentropic vortex patch, the
! amplitude of its disturbance and also its domain of influence
x_centroid = patch_icpp(patch_id)%x_centroid
y_centroid = patch_icpp(patch_id)%y_centroid
epsilon = patch_icpp(patch_id)%epsilon
beta = patch_icpp(patch_id)%beta

! Reference density, velocity, pressure and specific heat ratio
! function of the isentropic vortex patch
rho = patch_icpp(patch_id)%rho
vel = patch_icpp(patch_id)%vel
pres = patch_icpp(patch_id)%pres
gamma = patch_icpp(patch_id)%gamma

! Density
q_prim_vf(1)%sf(j, k, 0) = &
rho*(1d0 - (rho/pres)*(epsilon/(2d0*pi))* &
(epsilon/(8d0*beta*(gamma + 1d0)*pi))* &
exp(2d0*beta*(1d0 - (x_cc(j) - x_centroid)**2 &
- (y_cc(k) - y_centroid)**2)) &
)**gamma

! Velocity
q_prim_vf(2)%sf(j, k, 0) = &
vel(1) - (y_cc(k) - y_centroid)*(epsilon/(2d0*pi))* &
exp(beta*(1d0 - (x_cc(j) - x_centroid)**2 &
- (y_cc(k) - y_centroid)**2))
q_prim_vf(3)%sf(j, k, 0) = &
vel(2) + (x_cc(j) - x_centroid)*(epsilon/(2d0*pi))* &
exp(beta*(1d0 - (x_cc(j) - x_centroid)**2 &
- (y_cc(k) - y_centroid)**2))

! Pressure
q_prim_vf(4)%sf(j, k, 0) = &
pres*(1d0 - (rho/pres)*(epsilon/(2d0*pi))* &
(epsilon/(8d0*beta*(gamma + 1d0)*pi))* &
exp(2d0*beta*(1d0 - (x_cc(j) - x_centroid)**2 &
- (y_cc(k) - y_centroid)**2)) &
)**(gamma + 1d0)

! Specific heat ratio function
q_prim_vf(5)%sf(j, k, 0) = gamma

! Liquid stiffness function
q_prim_vf(6)%sf(j, k, 0) = 0d0
! Specific heat ratio function
q_prim_vf(gamma_idx)%sf(j, k, l) = &
eta*patch_icpp(patch_id)%gamma &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%gamma

end if
! Pressure
q_prim_vf(E_idx)%sf(j, k, l) = &
1d0/q_prim_vf(gamma_idx)%sf(j, k, l)* &
(eta*patch_icpp(patch_id)%gamma &
*patch_icpp(patch_id)%pres &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%gamma &
*patch_icpp(smooth_patch_id)%pres)

! Liquid stiffness function
q_prim_vf(pi_inf_idx)%sf(j, k, l) = &
eta*patch_icpp(patch_id)%pi_inf &
+ (1d0 - eta)*patch_icpp(smooth_patch_id)%pi_inf

! Updating the patch identities bookkeeping variable
if (1d0 - eta < 1d-16) patch_id_fp(j, k, l) = patch_id
Expand Down Expand Up @@ -598,47 +546,6 @@ subroutine s_assign_patch_species_primitive_variables(patch_id, j, k, l, &
end do
end if

if (patch_icpp(patch_id)%geometry == 6) then
x_centroid = patch_icpp(patch_id)%x_centroid
y_centroid = patch_icpp(patch_id)%y_centroid
epsilon = patch_icpp(patch_id)%epsilon
beta = patch_icpp(patch_id)%beta

! Reference density, velocity, pressure and specific heat ratio
! function of the isentropic vortex patch
rho = patch_icpp(patch_id)%rho
vel = patch_icpp(patch_id)%vel
pres = patch_icpp(patch_id)%pres
gamma = patch_icpp(patch_id)%gamma

! Density
q_prim_vf(1)%sf(j, k, 0) = &
rho*(1d0 - (rho/pres)*(epsilon/(2d0*pi))* &
(epsilon/(8d0*beta*(gamma + 1d0)*pi))* &
exp(2d0*beta*(1d0 - (x_cc(j) - x_centroid)**2 &
- (y_cc(k) - y_centroid)**2)) &
)**gamma

! Velocity
q_prim_vf(2)%sf(j, k, 0) = &
vel(1) - (y_cc(k) - y_centroid)*(epsilon/(2d0*pi))* &
exp(beta*(1d0 - (x_cc(j) - x_centroid)**2 &
- (y_cc(k) - y_centroid)**2))
q_prim_vf(3)%sf(j, k, 0) = &
vel(2) + (x_cc(j) - x_centroid)*(epsilon/(2d0*pi))* &
exp(beta*(1d0 - (x_cc(j) - x_centroid)**2 &
- (y_cc(k) - y_centroid)**2))

! Pressure
q_prim_vf(4)%sf(j, k, 0) = &
pres*(1d0 - (rho/pres)*(epsilon/(2d0*pi))* &
(epsilon/(8d0*beta*(gamma + 1d0)*pi))* &
exp(2d0*beta*(1d0 - (x_cc(j) - x_centroid)**2 &
- (y_cc(k) - y_centroid)**2)) &
)**(gamma + 1d0)

end if

! Updating the patch identities bookkeeping variable
if (1d0 - eta < 1d-16) patch_id_fp(j, k, l) = patch_id

Expand Down
30 changes: 3 additions & 27 deletions src/pre_process/m_check_patches.fpp
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,9 @@ contains
elseif (patch_icpp(i)%geometry == 5) then
call s_check_ellipse_patch_geometry(i)
elseif (patch_icpp(i)%geometry == 6) then
call s_check_isentropic_vortex_patch_geometry(i)
call s_mpi_abort('Unimplemented choice of geometry 6'// &
' (formerly "Vortex") of active patch '//trim(iStr)// &
' detected. Exiting ...')
elseif (patch_icpp(i)%geometry == 7) then
call s_check_2D_analytical_patch_geometry(i)
elseif (patch_icpp(i)%geometry == 8) then
Expand Down Expand Up @@ -270,32 +272,6 @@ contains

end subroutine s_check_ellipse_patch_geometry ! ------------------------

!> This subroutine verifies that the geometric parameters of
!! the isentropic vortex patch have been entered by the user
!! consistently.
!! @param patch_id Patch identifier
subroutine s_check_isentropic_vortex_patch_geometry(patch_id) ! --------

integer, intent(IN) :: patch_id
call s_int_to_str(patch_id, iStr)

! Constraints on the isentropic vortex patch geometric parameters
if (n == 0 .or. p > 0 .or. model_eqns == 2 &
.or. &
patch_icpp(patch_id)%radius <= 0d0 &
.or. &
patch_icpp(patch_id)%epsilon <= 0d0 &
.or. &
patch_icpp(patch_id)%beta <= 0d0) then

call s_mpi_abort('Inconsistency(ies) detected in '// &
'geometric parameters of isentropic '// &
'vortex patch '//trim(iStr)//'. Exiting ...')

end if

end subroutine s_check_isentropic_vortex_patch_geometry ! --------------

!> This subroutine verifies that the geometric parameters of
!! the Taylor Green vortex patch have been entered by the user
!! consistently.
Expand Down
Loading

0 comments on commit 61f84b3

Please sign in to comment.