diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 3c2fba720..d003e6fca 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -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`. | diff --git a/examples/2D_isentropicvortex/README.md b/examples/2D_isentropicvortex/README.md new file mode 100644 index 000000000..72032fbe1 --- /dev/null +++ b/examples/2D_isentropicvortex/README.md @@ -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) diff --git a/examples/2D_isentropicvortex/alpha_rho1.png b/examples/2D_isentropicvortex/alpha_rho1.png new file mode 100644 index 000000000..8793c2d02 Binary files /dev/null and b/examples/2D_isentropicvortex/alpha_rho1.png differ diff --git a/examples/2D_isentropicvortex/case.py b/examples/2D_isentropicvortex/case.py new file mode 100644 index 000000000..76606f4d9 --- /dev/null +++ b/examples/2D_isentropicvortex/case.py @@ -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, + # ========================================================================== +})) + +# ============================================================================== diff --git a/examples/2D_isentropicvortex/density_norms.png b/examples/2D_isentropicvortex/density_norms.png new file mode 100644 index 000000000..6cd066664 Binary files /dev/null and b/examples/2D_isentropicvortex/density_norms.png differ diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index f1bb1a7b9..4b634c3e9 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -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 diff --git a/src/pre_process/m_assign_variables.f90 b/src/pre_process/m_assign_variables.f90 index 74970aae3..db2f127c4 100644 --- a/src/pre_process/m_assign_variables.f90 +++ b/src/pre_process/m_assign_variables.f90 @@ -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 @@ -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 diff --git a/src/pre_process/m_check_patches.fpp b/src/pre_process/m_check_patches.fpp index 521e481a9..ec0421a6c 100644 --- a/src/pre_process/m_check_patches.fpp +++ b/src/pre_process/m_check_patches.fpp @@ -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 @@ -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. diff --git a/src/pre_process/m_initial_condition.fpp b/src/pre_process/m_initial_condition.fpp index 5fab877b7..c15be1c38 100644 --- a/src/pre_process/m_initial_condition.fpp +++ b/src/pre_process/m_initial_condition.fpp @@ -21,6 +21,10 @@ module m_initial_condition use m_global_parameters ! Global parameters for the code + use m_mpi_proxy !< Message passing interface (MPI) module proxy + + use m_helper + use m_variables_conversion ! Subroutines to change the state variables from ! one form to another @@ -28,8 +32,8 @@ module m_initial_condition use m_assign_variables - use m_eigen_solver ! Subroutines to solve eigenvalue problem for - ! complex general matrix + use m_eigen_solver ! Subroutines to solve eigenvalue problem for + ! complex general matrix ! ========================================================================== ! ========================================================================== @@ -107,6 +111,8 @@ contains integer :: i !< Generic loop operator + character(len=10) :: iStr + ! Converting the conservative variables to the primitive ones given ! preexisting initial condition data files were read in on start-up if (old_ic) then @@ -190,9 +196,10 @@ contains elseif (patch_icpp(i)%geometry == 5) then call s_ellipse(i, patch_id_fp, q_prim_vf) - ! Isentropic vortex patch + ! Unimplemented patch (formerly isentropic vortex) elseif (patch_icpp(i)%geometry == 6) then - call s_isentropic_vortex(i, patch_id_fp, q_prim_vf) + call s_mpi_abort('This used to be the isentropic vortex patch, '// & + 'which no longer exists. See Examples. Exiting ...') ! Analytical function patch for testing purposes elseif (patch_icpp(i)%geometry == 7) then diff --git a/src/pre_process/m_patches.fpp b/src/pre_process/m_patches.fpp index a561427c3..2f7329357 100644 --- a/src/pre_process/m_patches.fpp +++ b/src/pre_process/m_patches.fpp @@ -37,7 +37,6 @@ module m_patches s_ellipsoid, & s_rectangle, & s_sweep_line, & - s_isentropic_vortex, & s_2D_TaylorGreen_vortex, & s_1D_analytical, & s_1d_bubble_pulse, & @@ -646,61 +645,6 @@ contains end subroutine s_sweep_line ! ------------------------------------------ - !> The isentropic vortex is a 2D geometry that may be used, - !! for example, to generate an isentropic flow disturbance. - !! Geometry of the patch is well-defined when its centroid - !! and radius are provided. Notice that the patch DOES NOT - !! allow for the smoothing of its boundary. - !! @param patch_id is the patch identifier - subroutine s_isentropic_vortex(patch_id, patch_id_fp, q_prim_vf) ! ---------------------------- - - ! Patch identifier - integer, intent(IN) :: patch_id - integer, intent(INOUT), dimension(0:m, 0:n, 0:p) :: patch_id_fp - type(scalar_field), dimension(1:sys_size) :: q_prim_vf - - ! Generic loop iterators - integer :: i, j, k - - real(kind(0d0)) :: radius - - ! Transferring isentropic vortex patch's centroid and radius info - x_centroid = patch_icpp(patch_id)%x_centroid - y_centroid = patch_icpp(patch_id)%y_centroid - radius = patch_icpp(patch_id)%radius - - ! Since the isentropic vortex patch does not allow for its boundary - ! to get smoothed, the pseudo volume fraction is set to 1 to ensure - ! that only the current patch contributes to the fluid state in the - ! cells that this patch covers. - eta = 1d0 - - ! Verifying whether the isentropic vortex includes a particular cell - ! and verifying whether the current patch has permission to write to - ! that cell. If both queries work out the primitive variables of the - ! the current patch are assigned to this cell. - do j = 0, n - do i = 0, m - - if ((x_cc(i) - x_centroid)**2 & - + (y_cc(j) - y_centroid)**2 <= radius**2 & - .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, & - eta, q_prim_vf, patch_id_fp) - - @:analytical() - - end if - - end do - end do - - end subroutine s_isentropic_vortex ! ----------------------------------- - !> The Taylor Green vortex is 2D decaying vortex that may be used, !! for example, to verify the effects of viscous attenuation. !! Geometry of the patch is well-defined when its centroid