|
17 | 17 | #:enddef |
18 | 18 |
|
19 | 19 | #:def Hardcoded2D() |
20 | | - select case (patch_icpp(patch_id)%hcid) ! 2D_hardcoded_ic example case |
| 20 | + select case (patch_icpp(patch_id)%hcid) ! 2D_hardcoded_ic example case |
21 | 21 | case (200) |
22 | 22 | if (y_cc(j) <= (-x_cc(i)**3 + 1)**(1._wp/3._wp)) then |
23 | 23 | ! Volume Fractions |
24 | 24 | q_prim_vf(advxb)%sf(i, j, 0) = eps |
25 | 25 | q_prim_vf(advxe)%sf(i, j, 0) = 1._wp - eps |
26 | | - ! Densities |
27 | 26 | q_prim_vf(contxb)%sf(i, j, 0) = eps*1000._wp |
28 | 27 | q_prim_vf(contxe)%sf(i, j, 0) = (1._wp - eps)*1._wp |
29 | | - ! Pressure |
30 | 28 | q_prim_vf(E_idx)%sf(i, j, 0) = 1000._wp |
31 | 29 | end if |
32 | | - case (202) ! Gresho vortex (Gouasmi et al 2022 JCP) |
| 30 | + case (202) ! Gresho vortex (Gouasmi et al 2022 JCP) |
33 | 31 | r = ((x_cc(i) - 0.5_wp)**2 + (y_cc(j) - 0.5_wp)**2)**0.5_wp |
34 | 32 | rmax = 0.2_wp |
35 | 33 |
|
|
50 | 48 | q_prim_vf(momxe)%sf(i, j, 0) = 0._wp |
51 | 49 | q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2*(-2 + 4*log(2._wp)) |
52 | 50 | end if |
53 | | - case (203) ! Gresho vortex (Gouasmi et al 2022 JCP) with density correction |
| 51 | + case (203) ! Gresho vortex (Gouasmi et al 2022 JCP) with density correction |
54 | 52 | r = ((x_cc(i) - 0.5_wp)**2._wp + (y_cc(j) - 0.5_wp)**2)**0.5_wp |
55 | 53 | rmax = 0.2_wp |
56 | 54 |
|
|
73 | 71 | end if |
74 | 72 |
|
75 | 73 | q_prim_vf(contxb)%sf(i, j, 0) = q_prim_vf(E_idx)%sf(i, j, 0)**(1._wp/gam) |
76 | | - case (204) ! Rayleigh-Taylor instability |
| 74 | + case (204) ! Rayleigh-Taylor instability |
77 | 75 | rhoH = 3._wp |
78 | 76 | rhoL = 1._wp |
79 | 77 | pRef = 1.e5_wp |
|
104 | 102 | pInt = pref + rhoH*9.81_wp*(1.2_wp - intH) |
105 | 103 | q_prim_vf(E_idx)%sf(i, j, 0) = pInt + rhoL*9.81_wp*(intH - y_cc(j)) |
106 | 104 | end if |
107 | | - case (205) ! 2D lung wave interaction problem |
108 | | - h = 0.0_wp ! non dim origin y |
109 | | - lam = 1.0_wp ! non dim lambda |
110 | | - amp = patch_icpp(patch_id)%a(2) ! to be changed later! !non dim amplitude |
| 105 | + case (205) ! 2D lung wave interaction problem |
| 106 | + h = 0.0_wp ! non dim origin y |
| 107 | + lam = 1.0_wp ! non dim lambda |
| 108 | + amp = patch_icpp(patch_id)%a(2) ! to be changed later! !non dim amplitude |
111 | 109 |
|
112 | 110 | intH = amp*sin(2*pi*x_cc(i)/lam - pi/2) + h |
113 | 111 |
|
|
118 | 116 | q_prim_vf(advxb)%sf(i, j, 0) = patch_icpp(1)%alpha(1) |
119 | 117 | q_prim_vf(advxe)%sf(i, j, 0) = patch_icpp(1)%alpha(2) |
120 | 118 | end if |
121 | | - case (206) ! 2D lung wave interaction problem - horizontal domain |
122 | | - h = 0.0_wp ! non dim origin y |
123 | | - lam = 1.0_wp ! non dim lambda |
| 119 | + case (206) ! 2D lung wave interaction problem - horizontal domain |
| 120 | + h = 0.0_wp ! non dim origin y |
| 121 | + lam = 1.0_wp ! non dim lambda |
124 | 122 | amp = patch_icpp(patch_id)%a(2) |
125 | 123 |
|
126 | 124 | intL = amp*sin(2*pi*y_cc(j)/lam - pi/2) + h |
127 | 125 |
|
128 | | - if (x_cc(i) > intL) then ! this is the liquid |
| 126 | + if (x_cc(i) > intL) then ! this is the liquid |
129 | 127 | q_prim_vf(contxb)%sf(i, j, 0) = patch_icpp(1)%alpha_rho(1) |
130 | 128 | q_prim_vf(contxe)%sf(i, j, 0) = patch_icpp(1)%alpha_rho(2) |
131 | 129 | q_prim_vf(E_idx)%sf(i, j, 0) = patch_icpp(1)%pres |
132 | 130 | q_prim_vf(advxb)%sf(i, j, 0) = patch_icpp(1)%alpha(1) |
133 | 131 | q_prim_vf(advxe)%sf(i, j, 0) = patch_icpp(1)%alpha(2) |
134 | 132 | end if |
135 | | - case (207) ! Kelvin Helmholtz Instability |
| 133 | + case (207) ! Kelvin Helmholtz Instability |
136 | 134 | sigma = 0.05_wp/sqrt(2.0_wp) |
137 | 135 | gauss1 = exp(-(y_cc(j) - 0.75_wp)**2/(2.0_wp*sigma**2)) |
138 | 136 | gauss2 = exp(-(y_cc(j) - 0.25_wp)**2/(2.0_wp*sigma**2)) |
139 | 137 | q_prim_vf(momxb + 1)%sf(i, j, 0) = 0.1_wp*sin(4.0_wp*pi*x_cc(i))*(gauss1 + gauss2) |
140 | | - case (208) ! Richtmeyer Meshkov Instability |
| 138 | + case (208) ! Richtmeyer Meshkov Instability |
141 | 139 | lam = 1.0_wp |
142 | 140 | eps = 1.0e-6_wp |
143 | 141 | ei = 5.0_wp |
|
152 | 150 | q_prim_vf(advxb)%sf(i, j, 0) = alpha_sf6 |
153 | 151 | q_prim_vf(advxe)%sf(i, j, 0) = alpha_air |
154 | 152 | end if |
155 | | - case (250) ! MHD Orszag-Tang vortex |
| 153 | + case (250) ! MHD Orszag-Tang vortex |
156 | 154 | ! gamma = 5/3 rho = 25/(36*pi) p = 5/(12*pi) v = (-sin(2*pi*y), sin(2*pi*x), 0) B = (-sin(2*pi*y)/sqrt(4*pi), |
157 | 155 | ! sin(4*pi*x)/sqrt(4*pi), 0) |
158 | 156 |
|
|
161 | 159 |
|
162 | 160 | q_prim_vf(B_idx%beg)%sf(i, j, 0) = -sin(2._wp*pi*y_cc(j))/sqrt(4._wp*pi) |
163 | 161 | q_prim_vf(B_idx%beg + 1)%sf(i, j, 0) = sin(4._wp*pi*x_cc(i))/sqrt(4._wp*pi) |
164 | | - case (251) ! RMHD Cylindrical Blast Wave [Mignone, 2006: Section 4.3.1] |
| 162 | + case (251) ! RMHD Cylindrical Blast Wave [Mignone, 2006: Section 4.3.1] |
165 | 163 | if (x_cc(i)**2 + y_cc(j)**2 < 0.08_wp**2) then |
166 | 164 | q_prim_vf(contxb)%sf(i, j, 0) = 0.01 |
167 | 165 | q_prim_vf(E_idx)%sf(i, j, 0) = 1.0 |
|
176 | 174 | end if |
177 | 175 |
|
178 | 176 | ! case 252 is for the 2D MHD Rotor problem |
179 | | - case (252) ! 2D MHD Rotor Problem |
| 177 | + case (252) ! 2D MHD Rotor Problem |
180 | 178 | ! Ambient conditions are set in the JSON file. This case imposes the dense, rotating cylinder. |
181 | 179 | ! |
182 | 180 | ! gamma = 1.4 Ambient medium (r > 0.1): rho = 1, p = 1, v = 0, B = (1,0,0) Rotor (r <= 0.1): rho = 10, p = 1 v has angular |
|
202 | 200 | q_prim_vf(momxb)%sf(i, j, 0) = -(2._wp/sqrt(r_sq))*(y_cc(j) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp) |
203 | 201 | q_prim_vf(momxb + 1)%sf(i, j, 0) = (2._wp/sqrt(r_sq))*(x_cc(i) - 0.5_wp)*(0.115_wp - sqrt(r_sq))/(0.015_wp) |
204 | 202 | end if |
205 | | - case (253) ! MHD Smooth Magnetic Vortex |
| 203 | + case (253) ! MHD Smooth Magnetic Vortex |
206 | 204 | ! Section 5.2 of Implicit hybridized discontinuous Galerkin methods for compressible magnetohydrodynamics C. Ciuca, P. |
207 | 205 | ! Fernandez, A. Christophe, N.C. Nguyen, J. Peraire |
208 | 206 |
|
|
217 | 215 | ! pressure |
218 | 216 | q_prim_vf(E_idx)%sf(i, j, & |
219 | 217 | & 0) = 1._wp + (1 - 2._wp*(x_cc(i)**2 + y_cc(j)**2))*exp(1 - (x_cc(i)**2 + y_cc(j)**2))/((2._wp*pi)**3) |
220 | | - case (260) ! Gaussian Divergence Pulse |
| 218 | + case (260) ! Gaussian Divergence Pulse |
221 | 219 | ! Bx(x) = 1 + C * erf((x-0.5)/\sigma) => \partialBx/\partialx = C * (2/\sqrt\pi) * exp[-((x-0.5)/\sigma)**2] * (1/\sigma) |
222 | 220 | ! Choose C = \epsilon * \sigma * \sqrt\pi / 2 => \partialBx/\partialx = \epsilon * exp[-((x-0.5)/\sigma)**2] \psi is |
223 | 221 | ! initialized to zero everywhere. |
|
228 | 226 |
|
229 | 227 | ! B-field |
230 | 228 | q_prim_vf(B_idx%beg)%sf(i, j, 0) = 1._wp + C_mhd*erf((x_cc(i) - 0.5_wp)/sigma) |
231 | | - case (261) ! Blob |
| 229 | + case (261) ! Blob |
232 | 230 | r0 = 1._wp/sqrt(8._wp) |
233 | 231 | r2 = x_cc(i)**2 + y_cc(j)**2 |
234 | 232 | r = sqrt(r2) |
|
239 | 237 | ! q_prim_vf(B_idx%beg)%sf(i,j,0) = 1._wp/(4._wp*pi) * (alpha**8 - 2._wp*alpha**4 + 1._wp) q_prim_vf(E_idx)%sf(i,j,0) = |
240 | 238 | ! 6._wp - q_prim_vf(B_idx%beg)%sf(i,j,0)**2/2._wp |
241 | 239 | end if |
242 | | - case (262) ! Tilted 2D MHD shock‐tube at α = arctan2 (≈63.4°) |
| 240 | + case (262) ! Tilted 2D MHD shock‐tube at α = arctan2 (≈63.4°) |
243 | 241 | ! rotate by \alpha = atan(2) |
244 | 242 | alpha = atan(2._wp) |
245 | 243 | cosA = cos(alpha) |
|
0 commit comments