Skip to content

Commit 4ae6372

Browse files
sbryngelsonsbryngelson
authored
Merge branch 'master' into IGRMergeOne
2 parents 545c7c7 + da8ae1c commit 4ae6372

File tree

3 files changed

+68
-127
lines changed

3 files changed

+68
-127
lines changed

‎benchmarks/ibm/case.py‎

Lines changed: 32 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@
2323

2424
dx = 1.0 / (1.0 * (Nx + 1))
2525

26-
Tend = 64e-06
26+
Tend = 1e-4
2727
Nt = 200
2828
mydt = Tend / (1.0 * Nt)
2929

@@ -48,21 +48,21 @@
4848
"t_step_stop": int(20 * (95 * size + 5)),
4949
"t_step_save": int(20 * (95 * size + 5)),
5050
# Simulation Algorithm Parameters
51-
"num_patches": 2,
51+
"num_patches": 1,
5252
"model_eqns": 2,
5353
"alt_soundspeed": "F",
54-
"num_fluids": 1,
54+
"num_fluids": 2,
5555
"mpp_lim": "F",
5656
"mixture_err": "F",
5757
"time_stepper": 3,
5858
"weno_order": 5,
5959
"weno_eps": 1.0e-16,
60-
"weno_Re_flux": "F",
61-
"weno_avg": "F",
62-
"mapped_weno": "F",
60+
"weno_Re_flux": "T",
61+
"weno_avg": "T",
62+
"mapped_weno": "T",
6363
"null_weights": "F",
6464
"mp_weno": "F",
65-
"riemann_solver": 1,
65+
"riemann_solver": 2,
6666
"wave_speeds": 1,
6767
"avg_state": 2,
6868
"bc_x%beg": -3,
@@ -71,48 +71,45 @@
7171
"bc_y%end": -3,
7272
"bc_z%beg": -3,
7373
"bc_z%end": -3,
74-
# Turning on Hypoelasticity
75-
"hypoelasticity": "T",
76-
"fd_order": 4,
74+
# Turn on IBM
75+
"ib": "T",
76+
"num_ibs": 1,
77+
"viscous": "T",
7778
# Formatted Database Files Structure Parameters
7879
"format": 1,
7980
"precision": 2,
8081
"prim_vars_wrt": "T",
8182
"parallel_io": "T",
8283
# Patch 1 L
8384
"patch_icpp(1)%geometry": 9,
84-
"patch_icpp(1)%x_centroid": 0.25,
85+
"patch_icpp(1)%x_centroid": 0.5,
8586
"patch_icpp(1)%y_centroid": 0.25,
8687
"patch_icpp(1)%z_centroid": 0.25,
87-
"patch_icpp(1)%length_x": 0.5,
88+
"patch_icpp(1)%length_x": 1.0,
8889
"patch_icpp(1)%length_y": 0.5,
8990
"patch_icpp(1)%length_z": 0.5,
90-
"patch_icpp(1)%vel(1)": 0.0,
91+
"patch_icpp(1)%vel(1)": 0.1,
9192
"patch_icpp(1)%vel(2)": 0,
9293
"patch_icpp(1)%vel(3)": 0,
93-
"patch_icpp(1)%pres": 1.0e8,
94-
"patch_icpp(1)%alpha_rho(1)": 1000,
95-
"patch_icpp(1)%alpha(1)": 1.0,
96-
"patch_icpp(1)%tau_e(1)": 0.0,
97-
# Patch 2 R
98-
"patch_icpp(2)%geometry": 9,
99-
"patch_icpp(2)%x_centroid": 0.75,
100-
"patch_icpp(2)%y_centroid": 0.25,
101-
"patch_icpp(2)%z_centroid": 0.25,
102-
"patch_icpp(2)%length_x": 0.5,
103-
"patch_icpp(2)%length_y": 0.5,
104-
"patch_icpp(2)%length_z": 0.5,
105-
"patch_icpp(2)%vel(1)": 0,
106-
"patch_icpp(2)%vel(2)": 0,
107-
"patch_icpp(2)%vel(3)": 0,
108-
"patch_icpp(2)%pres": 1.0e05,
109-
"patch_icpp(2)%alpha_rho(1)": 1000,
110-
"patch_icpp(2)%alpha(1)": 1.0,
111-
"patch_icpp(2)%tau_e(1)": 0.0,
94+
"patch_icpp(1)%pres": 1.0,
95+
"patch_icpp(1)%alpha_rho(1)": 0.8e00,
96+
"patch_icpp(1)%alpha(1)": 0.8e00,
97+
"patch_icpp(1)%alpha_rho(2)": 0.2e00,
98+
"patch_icpp(1)%alpha(2)": 0.2e00,
99+
# Patch: Sphere Immersed Boundary
100+
"patch_ib(1)%geometry": 8,
101+
"patch_ib(1)%x_centroid": 0.25,
102+
"patch_ib(1)%y_centroid": 0.25,
103+
"patch_ib(1)%z_centroid": 0.25,
104+
"patch_ib(1)%radius": 0.1,
112105
# Fluids Physical Parameters
113-
"fluid_pp(1)%gamma": 1.0e00 / (4.4e00 - 1.0e00),
114-
"fluid_pp(1)%pi_inf": 4.4e00 * 6.0e08 / (4.4e00 - 1.0e00),
115-
"fluid_pp(1)%G": 10e09,
106+
# Specify 2 fluids
107+
"fluid_pp(1)%gamma": 1.0e00 / (1.4 - 1.0e00),
108+
"fluid_pp(1)%pi_inf": 0,
109+
"fluid_pp(1)%Re(1)": 54000,
110+
"fluid_pp(2)%gamma": 1.0e00 / (1.4 - 1.0e00),
111+
"fluid_pp(2)%pi_inf": 0,
112+
"fluid_pp(2)%Re(1)": 54000,
116113
}
117114
)
118115
)

‎src/post_process/m_mpi_proxy.fpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -108,6 +108,7 @@ contains
108108
& 'adv_n', 'ib', 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', &
109109
& 'surface_tension', 'hyperelasticity', 'bubbles_lagrange', &
110110
& 'output_partial_domain', 'relativity', 'cont_damage', 'bc_io' ]
111+
111112
call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
112113
#:endfor
113114

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 35 additions & 92 deletions
Original file line numberDiff line numberDiff line change
@@ -2995,19 +2995,15 @@ contains
29952995
! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
29962996
if (mhd) then
29972997
if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
2998-
B%L(1) = Bx0
2999-
B%R(1) = Bx0
3000-
B%L(2) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg)
3001-
B%R(2) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg)
3002-
B%L(3) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1)
3003-
B%R(3) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + 1)
2998+
B%L = [Bx0, qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg), qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1)]
2999+
B%R = [Bx0, qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg), qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + 1)]
30043000
else ! 2D/3D: Bx, By, Bz as variables
3005-
B%L(1) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1) - 1)
3006-
B%R(1) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(1) - 1)
3007-
B%L(2) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1)
3008-
B%R(2) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(2) - 1)
3009-
B%L(3) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1)
3010-
B%R(3) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(3) - 1)
3001+
B%L = [qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1) - 1), &
3002+
qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1), &
3003+
qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1)]
3004+
B%R = [qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(1) - 1), &
3005+
qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(2) - 1), &
3006+
qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(3) - 1)]
30113007
end if
30123008
end if
30133009

@@ -3058,74 +3054,38 @@ contains
30583054
E_starL = ((s_L - vel%L(1))*E%L - pTot_L*vel%L(1) + p_star*s_M)/(s_L - s_M)
30593055
E_starR = ((s_R - vel%R(1))*E%R - pTot_R*vel%R(1) + p_star*s_M)/(s_R - s_M)
30603056

3061-
! (5) Compute the left/right conserved state vectors
3062-
U_L(1) = rho%L
3063-
U_L(2) = rho%L*vel%L(1)
3064-
U_L(3) = rho%L*vel%L(2)
3065-
U_L(4) = rho%L*vel%L(3)
3066-
U_L(5) = B%L(2)
3067-
U_L(6) = B%L(3)
3068-
U_L(7) = E%L
3069-
3070-
U_R(1) = rho%R
3071-
U_R(2) = rho%R*vel%R(1)
3072-
U_R(3) = rho%R*vel%R(2)
3073-
U_R(4) = rho%R*vel%R(3)
3074-
U_R(5) = B%R(2)
3075-
U_R(6) = B%R(3)
3076-
U_R(7) = E%R
3077-
3078-
! (6) Compute the left/right star state vectors
3079-
U_starL(1) = rhoL_star
3080-
U_starL(2) = rhoL_star*s_M
3081-
U_starL(3) = rhoL_star*vel%L(2)
3082-
U_starL(4) = rhoL_star*vel%L(3)
3083-
U_starL(5) = B%L(2)
3084-
U_starL(6) = B%L(3)
3085-
U_starL(7) = E_starL
3086-
3087-
U_starR(1) = rhoR_star
3088-
U_starR(2) = rhoR_star*s_M
3089-
U_starR(3) = rhoR_star*vel%R(2)
3090-
U_starR(4) = rhoR_star*vel%R(3)
3091-
U_starR(5) = B%R(2)
3092-
U_starR(6) = B%R(3)
3093-
U_starR(7) = E_starR
3094-
3095-
! (7) Compute the left/right fluxes
3096-
F_L(1) = rho%L*vel%L(1)
3097-
F_L(2) = rho%L*vel%L(1)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
3098-
F_L(3) = rho%L*vel%L(1)*vel%L(2) - B%L(1)*B%L(2)
3099-
F_L(4) = rho%L*vel%L(1)*vel%L(3) - B%L(1)*B%L(3)
3100-
F_L(5) = vel%L(1)*B%L(2) - vel%L(2)*B%L(1)
3101-
F_L(6) = vel%L(1)*B%L(3) - vel%L(3)*B%L(1)
3057+
! (5) Compute left/right state vectors and fluxes
3058+
U_L = [rho%L, rho%L*vel%L(1:3), B%L(2:3), E%L]
3059+
U_starL = [rhoL_star, rhoL_star*s_M, rhoL_star*vel%L(2:3), B%L(2:3), E_starL]
3060+
U_R = [rho%R, rho%R*vel%R(1:3), B%R(2:3), E%R]
3061+
U_starR = [rhoR_star, rhoR_star*s_M, rhoR_star*vel%R(2:3), B%R(2:3), E_starR]
3062+
3063+
! Compute the left/right fluxes
3064+
F_L(1) = U_L(2)
3065+
F_L(2) = U_L(2)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
3066+
F_L(3:4) = U_L(2)*vel%L(2:3) - B%L(1)*B%L(2:3)
3067+
F_L(5:6) = vel%L(1)*B%L(2:3) - vel%L(2:3)*B%L(1)
31023068
F_L(7) = (E%L + pTot_L)*vel%L(1) - B%L(1)*(vel%L(1)*B%L(1) + vel%L(2)*B%L(2) + vel%L(3)*B%L(3))
31033069

3104-
F_R(1) = rho%R*vel%R(1)
3105-
F_R(2) = rho%R*vel%R(1)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
3106-
F_R(3) = rho%R*vel%R(1)*vel%R(2) - B%R(1)*B%R(2)
3107-
F_R(4) = rho%R*vel%R(1)*vel%R(3) - B%R(1)*B%R(3)
3108-
F_R(5) = vel%R(1)*B%R(2) - vel%R(2)*B%R(1)
3109-
F_R(6) = vel%R(1)*B%R(3) - vel%R(3)*B%R(1)
3070+
F_R(1) = U_R(2)
3071+
F_R(2) = U_R(2)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
3072+
F_R(3:4) = U_R(2)*vel%R(2:3) - B%R(1)*B%R(2:3)
3073+
F_R(5:6) = vel%R(1)*B%R(2:3) - vel%R(2:3)*B%R(1)
31103074
F_R(7) = (E%R + pTot_R)*vel%R(1) - B%R(1)*(vel%R(1)*B%R(1) + vel%R(2)*B%R(2) + vel%R(3)*B%R(3))
3111-
3112-
! (8) Compute the left/right star fluxes (note array operations)
3075+
! Compute the star flux using HLL relation
31133076
F_starL = F_L + s_L*(U_starL - U_L)
31143077
F_starR = F_R + s_R*(U_starR - U_R)
3115-
3116-
! (9) Compute the rotational (Alfvén) speeds
3078+
! Compute the rotational (Alfvén) speeds
31173079
s_starL = s_M - abs(B%L(1))/sqrt(rhoL_star)
31183080
s_starR = s_M + abs(B%L(1))/sqrt(rhoR_star)
3081+
! Compute the double–star states [Miyoshi Eqns. (59)-(62)]
3082+
sqrt_rhoL_star = sqrt(rhoL_star); sqrt_rhoR_star = sqrt(rhoR_star)
3083+
vL_star = vel%L(2); wL_star = vel%L(3)
3084+
vR_star = vel%R(2); wR_star = vel%R(3)
31193085

3120-
! (10) Compute the double–star states [Miyoshi Eqns. (59)-(62)]
3121-
sqrt_rhoL_star = sqrt(rhoL_star)
3122-
sqrt_rhoR_star = sqrt(rhoR_star)
3086+
! (6) Compute the double–star states [Miyoshi Eqns. (59)-(62)]
31233087
denom_ds = sqrt_rhoL_star + sqrt_rhoR_star
31243088
sign_Bx = sign(1._wp, B%L(1))
3125-
vL_star = vel%L(2)
3126-
wL_star = vel%L(3)
3127-
vR_star = vel%R(2)
3128-
wR_star = vel%R(3)
31293089
v_double = (sqrt_rhoL_star*vL_star + sqrt_rhoR_star*vR_star + (B%R(2) - B%L(2))*sign_Bx)/denom_ds
31303090
w_double = (sqrt_rhoL_star*wL_star + sqrt_rhoR_star*wR_star + (B%R(3) - B%L(3))*sign_Bx)/denom_ds
31313091
By_double = (sqrt_rhoL_star*B%R(2) + sqrt_rhoR_star*B%L(2) + sqrt_rhoL_star*sqrt_rhoR_star*(vR_star - vL_star)*sign_Bx)/denom_ds
@@ -3135,21 +3095,8 @@ contains
31353095
E_doubleR = E_starR + sqrt_rhoR_star*((vR_star*B%R(2) + wR_star*B%R(3)) - (v_double*By_double + w_double*Bz_double))*sign_Bx
31363096
E_double = 0.5_wp*(E_doubleL + E_doubleR)
31373097

3138-
U_doubleL(1) = rhoL_star
3139-
U_doubleL(2) = rhoL_star*s_M
3140-
U_doubleL(3) = rhoL_star*v_double
3141-
U_doubleL(4) = rhoL_star*w_double
3142-
U_doubleL(5) = By_double
3143-
U_doubleL(6) = Bz_double
3144-
U_doubleL(7) = E_double
3145-
3146-
U_doubleR(1) = rhoR_star
3147-
U_doubleR(2) = rhoR_star*s_M
3148-
U_doubleR(3) = rhoR_star*v_double
3149-
U_doubleR(4) = rhoR_star*w_double
3150-
U_doubleR(5) = By_double
3151-
U_doubleR(6) = Bz_double
3152-
U_doubleR(7) = E_double
3098+
U_doubleL = [rhoL_star, rhoL_star*s_M, rhoL_star*v_double, rhoL_star*w_double, By_double, Bz_double, E_double]
3099+
U_doubleR = [rhoR_star, rhoR_star*s_M, rhoR_star*v_double, rhoR_star*w_double, By_double, Bz_double, E_double]
31533100

31543101
! (11) Choose HLLD flux based on wave-speed regions
31553102
if (0.0_wp <= s_L) then
@@ -3170,16 +3117,12 @@ contains
31703117
! Mass
31713118
flux_rs${XYZ}$_vf(j, k, l, 1) = F_hlld(1) ! TODO multi-component
31723119
! Momentum
3173-
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(1)) = F_hlld(2)
3174-
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(2)) = F_hlld(3)
3175-
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(3)) = F_hlld(4)
3120+
flux_rs${XYZ}$_vf(j, k, l, [contxe + dir_idx(1), contxe + dir_idx(2), contxe + dir_idx(3)]) = F_hlld([2, 3, 4])
31763121
! Magnetic field
31773122
if (n == 0) then
3178-
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg) = F_hlld(5)
3179-
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1) = F_hlld(6)
3123+
flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg, B_idx%beg + 1]) = F_hlld([5, 6])
31803124
else
3181-
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1) = F_hlld(5)
3182-
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1) = F_hlld(6)
3125+
flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg + dir_idx(2) - 1, B_idx%beg + dir_idx(3) - 1]) = F_hlld([5, 6])
31833126
end if
31843127
! Energy
31853128
flux_rs${XYZ}$_vf(j, k, l, E_idx) = F_hlld(7)

0 commit comments

Comments
 (0)