|
25 | 25 | end if
|
26 | 26 | case (202) ! Gresho vortex (Gouasmi et al 2022 JCP)
|
27 | 27 | r = ((x_cc(i) - 0.5_wp)**2 + (y_cc(j) - 0.5_wp)**2)**0.5_wp
|
28 |
| - rmax = 0.2 |
| 28 | + rmax = 0.2_wp |
29 | 29 |
|
30 | 30 | gam = 1._wp + 1._wp/fluid_pp(1)%gamma
|
31 | 31 | umax = 2*pi*rmax*patch_icpp(patch_id)%vel(2)
|
|
42 | 42 | else
|
43 | 43 | q_prim_vf(momxb)%sf(i, j, 0) = 0._wp
|
44 | 44 | q_prim_vf(momxe)%sf(i, j, 0) = 0._wp
|
45 |
| - q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2*(-2 + 4*log(2.)) |
| 45 | + q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2*(-2 + 4*log(2._wp)) |
46 | 46 | end if
|
47 | 47 | case (203) ! Gresho vortex (Gouasmi et al 2022 JCP) with density correction
|
48 |
| - r = ((x_cc(i) - 0.5_wp)**2 + (y_cc(j) - 0.5_wp)**2)**0.5_wp |
49 |
| - rmax = 0.2 |
| 48 | + r = ((x_cc(i) - 0.5_wp)**2._wp + (y_cc(j) - 0.5_wp)**2)**0.5_wp |
| 49 | + rmax = 0.2_wp |
50 | 50 |
|
51 | 51 | gam = 1._wp + 1._wp/fluid_pp(1)%gamma
|
52 | 52 | umax = 2*pi*rmax*patch_icpp(patch_id)%vel(2)
|
|
55 | 55 | if (r < rmax) then
|
56 | 56 | q_prim_vf(momxb)%sf(i, j, 0) = -(y_cc(j) - 0.5_wp)*umax/rmax
|
57 | 57 | q_prim_vf(momxe)%sf(i, j, 0) = (x_cc(i) - 0.5_wp)*umax/rmax
|
58 |
| - q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2*((r/rmax)**2/2._wp) |
| 58 | + q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2*((r/rmax)**2._wp/2._wp) |
59 | 59 | else if (r < 2*rmax) then
|
60 | 60 | q_prim_vf(momxb)%sf(i, j, 0) = -((y_cc(j) - 0.5_wp)/r)*umax*(2._wp - r/rmax)
|
61 | 61 | q_prim_vf(momxe)%sf(i, j, 0) = ((x_cc(i) - 0.5_wp)/r)*umax*(2._wp - r/rmax)
|
62 |
| - q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2*((r/rmax)**2/2._wp + 4*(1 - (r/rmax) + log(r/rmax))) |
| 62 | + q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2*((r/rmax)**2/2._wp + 4._wp*(1._wp - (r/rmax) + log(r/rmax))) |
63 | 63 | else
|
64 | 64 | q_prim_vf(momxb)%sf(i, j, 0) = 0._wp
|
65 | 65 | q_prim_vf(momxe)%sf(i, j, 0) = 0._wp
|
66 |
| - q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2*(-2 + 4*log(2.)) |
| 66 | + q_prim_vf(E_idx)%sf(i, j, 0) = p0 + umax**2._wp*(-2._wp + 4*log(2._wp)) |
67 | 67 | end if
|
68 | 68 |
|
69 | 69 | q_prim_vf(contxb)%sf(i, j, 0) = q_prim_vf(E_idx)%sf(i, j, 0)**(1._wp/gam)
|
70 | 70 |
|
71 | 71 | case (204) ! Rayleigh-Taylor instability
|
72 |
| - rhoH = 3 |
73 |
| - rhoL = 1 |
74 |
| - pRef = 1e5_wp |
| 72 | + rhoH = 3._wp |
| 73 | + rhoL = 1._wp |
| 74 | + pRef = 1.e5_wp |
75 | 75 | pInt = pRef
|
76 |
| - h = 0.7 |
77 |
| - lam = 0.2 |
78 |
| - wl = 2*pi/lam |
79 |
| - amp = 0.05/wl |
| 76 | + h = 0.7_wp |
| 77 | + lam = 0.2_wp |
| 78 | + wl = 2._wp*pi/lam |
| 79 | + amp = 0.05_wp/wl |
80 | 80 |
|
81 |
| - intH = amp*sin(2*pi*x_cc(i)/lam - pi/2) + h |
| 81 | + intH = amp*sin(2._wp*pi*x_cc(i)/lam - pi/2._wp) + h |
82 | 82 |
|
83 |
| - alph = 5e-1_wp*(1 + tanh((y_cc(j) - intH)/2.5e-3_wp)) |
| 83 | + alph = 0.5_wp*(1._wp + tanh((y_cc(j) - intH)/2.5e-3_wp)) |
84 | 84 |
|
85 | 85 | if (alph < eps) alph = eps
|
86 |
| - if (alph > 1 - eps) alph = 1 - eps |
| 86 | + if (alph > 1._wp - eps) alph = 1._wp - eps |
87 | 87 |
|
88 | 88 | if (y_cc(j) > intH) then
|
89 | 89 | q_prim_vf(advxb)%sf(i, j, 0) = alph
|
90 |
| - q_prim_vf(advxe)%sf(i, j, 0) = 1 - alph |
| 90 | + q_prim_vf(advxe)%sf(i, j, 0) = 1._wp - alph |
91 | 91 | q_prim_vf(contxb)%sf(i, j, 0) = alph*rhoH
|
92 |
| - q_prim_vf(contxe)%sf(i, j, 0) = (1 - alph)*rhoL |
93 |
| - q_prim_vf(E_idx)%sf(i, j, 0) = pref + rhoH*9.81*(1.2 - y_cc(j)) |
| 92 | + q_prim_vf(contxe)%sf(i, j, 0) = (1._wp - alph)*rhoL |
| 93 | + q_prim_vf(E_idx)%sf(i, j, 0) = pref + rhoH*9.81_wp*(1.2_wp - y_cc(j)) |
94 | 94 | else
|
95 | 95 | q_prim_vf(advxb)%sf(i, j, 0) = alph
|
96 |
| - q_prim_vf(advxe)%sf(i, j, 0) = 1 - alph |
| 96 | + q_prim_vf(advxe)%sf(i, j, 0) = 1._wp - alph |
97 | 97 | q_prim_vf(contxb)%sf(i, j, 0) = alph*rhoH
|
98 |
| - q_prim_vf(contxe)%sf(i, j, 0) = (1 - alph)*rhoL |
99 |
| - pInt = pref + rhoH*9.81*(1.2 - intH) |
100 |
| - q_prim_vf(E_idx)%sf(i, j, 0) = pInt + rhoL*9.81*(intH - y_cc(j)) |
| 98 | + q_prim_vf(contxe)%sf(i, j, 0) = (1._wp - alph)*rhoL |
| 99 | + pInt = pref + rhoH*9.81_wp*(1.2_wp - intH) |
| 100 | + q_prim_vf(E_idx)%sf(i, j, 0) = pInt + rhoL*9.81_wp*(intH - y_cc(j)) |
101 | 101 | end if
|
102 | 102 |
|
103 | 103 | case default
|
|
0 commit comments