Repository navigation
2nd order IBM #1958
+657
−55
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,119 @@ | ||
| #!/usr/bin/env python3 | ||
| """ | ||
| NACA 0012 at M = 0.3, alpha = 2 deg, 200 cells per chord: the 2D_ibm_airfoil_surface_pressure example | ||
| switched to a viscous no-slip wall (Re = 1000 per chord), with frequent snapshots for time averaging. | ||
| Env: RE, T_STOP, T_SAVE, SECOND_ORDER=1 (PR branch only: second-order IB velocities). | ||
| """ | ||
|
|
||
| import json | ||
| import math | ||
| import os | ||
|
|
||
| Ma = 0.3 | ||
| alpha_deg = 2.0 | ||
| gamma = 1.4 | ||
| rho_inf, U_inf, chord = 1.0, 1.0, 1.0 | ||
| Re = float(os.environ.get("RE", 1000.0)) | ||
| P_inf = rho_inf * U_inf**2 / (gamma * Ma**2) | ||
|
|
||
| case = { | ||
| # --- Output --- | ||
| "run_time_info": "T", | ||
| "format": 2, | ||
| "precision": 2, | ||
| "parallel_io": "T", | ||
| "prim_vars_wrt": "T", | ||
| # --- Domain: pre-stretch extents; the stretching maps them outward --- | ||
| "x_domain%beg": -3.0, | ||
| "x_domain%end": 4.0, | ||
| "y_domain%beg": -3.0, | ||
| "y_domain%end": 3.0, | ||
| "m": 1399, | ||
| "n": 1199, | ||
| "p": 0, | ||
| "cyl_coord": "F", | ||
| "stretch_x": "T", | ||
| "a_x": 15.0, | ||
| "x_a": -0.8, | ||
| "x_b": 1.8, | ||
| "loops_x": 2, | ||
| "stretch_y": "T", | ||
| "a_y": 15.0, | ||
| "y_a": -0.7, | ||
| "y_b": 0.7, | ||
| "loops_y": 2, | ||
| # --- Time stepping --- | ||
| "cfl_adap_dt": "T", | ||
| "cfl_target": 0.5, | ||
| "n_start": 0, | ||
| "t_save": 0.2, | ||
| "t_stop": 6.0, | ||
| # --- Numerics --- | ||
| "num_patches": 1, | ||
| "num_fluids": 1, | ||
| "model_eqns": 2, | ||
| "alt_soundspeed": "F", | ||
| "mpp_lim": "F", | ||
| "mixture_err": "T", | ||
| "time_stepper": 3, | ||
| "weno_order": 5, | ||
| "weno_eps": 1.0e-10, | ||
| "weno_Re_flux": "F", | ||
| "weno_avg": "T", | ||
| "avg_state": 2, | ||
| "mapped_weno": "T", | ||
| "null_weights": "F", | ||
| "mp_weno": "F", | ||
| "riemann_solver": 2, | ||
| "low_Mach": 2, | ||
| "wave_speeds": 1, | ||
| "viscous": "T", | ||
| "fd_order": 4, | ||
| # --- Uniform freestream --- | ||
| "patch_icpp(1)%geometry": 3, | ||
| "patch_icpp(1)%x_centroid": 0.0, | ||
| "patch_icpp(1)%y_centroid": 0.0, | ||
| "patch_icpp(1)%length_x": 1.0e3, | ||
| "patch_icpp(1)%length_y": 1.0e3, | ||
| "patch_icpp(1)%vel(1)": U_inf, | ||
| "patch_icpp(1)%vel(2)": 0.0, | ||
| "patch_icpp(1)%pres": P_inf, | ||
| "patch_icpp(1)%alpha_rho(1)": rho_inf, | ||
| "patch_icpp(1)%alpha(1)": 1.0, | ||
| "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0), | ||
| "fluid_pp(1)%eos": "ideal_gas", | ||
| "fluid_pp(1)%Re(1)": Re, | ||
| # --- Characteristic far-field boundaries --- | ||
| "bc_x%beg": -7, | ||
| "bc_x%grcbc_in": "T", | ||
| "bc_x%vel_in(1)": U_inf, | ||
| "bc_x%vel_in(2)": 0.0, | ||
| "bc_x%pres_in": P_inf, | ||
| "bc_x%alpha_rho_in(1)": rho_inf, | ||
| "bc_x%alpha_in(1)": 1.0, | ||
| "bc_x%end": -8, | ||
| "bc_x%grcbc_out": "T", | ||
| "bc_x%pres_out": P_inf, | ||
| "bc_y%beg": -9, | ||
| "bc_y%end": -9, | ||
| # --- IB: NACA 0012 (m must be > 0, so a negligible camber stands in for 0) --- | ||
| "ib": "T", | ||
| "num_ibs": 1, | ||
| "patch_ib(1)%geometry": 4, | ||
| "patch_ib(1)%x_centroid": 0.0, | ||
| "patch_ib(1)%y_centroid": 0.0, | ||
| "patch_ib(1)%airfoil_id": 1, | ||
| "patch_ib(1)%angles(3)": -math.radians(alpha_deg), | ||
| "patch_ib(1)%slip": "F", | ||
| "patch_ib(1)%moving_ibm": 0, | ||
| "ib_airfoil(1)%c": chord, | ||
| "ib_airfoil(1)%t": 0.12, | ||
| "ib_airfoil(1)%p": 0.4, | ||
| "ib_airfoil(1)%m": 1.0e-9, | ||
| } | ||
|
|
||
| # the parameters for second-order IBM | ||
| case.update({"ib_second_order_vel": "T", "ib_ip_min_dist": 1.5}) | ||
|
|
||
|
|
||
| print(json.dumps(case, indent=4)) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -30,7 +30,7 @@ module m_ibm | |
| implicit none | ||
|
|
||
| private :: s_compute_image_points, s_compute_interpolation_coeffs, s_interpolate_image_point, s_find_ghost_points, & | ||
| & s_find_num_ghost_points, s_compute_ghost_point_pressure, s_compute_ghost_point_velocity, s_interpolation_weights | ||
| ; public :: s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module, s_report_ibm_surface | ||
|
|
||
| !> Ghost points at which the reacting-surface Newton solve did not reach its tolerance, so the point fell back to a chemically | ||
|
|
@@ -83,6 +83,7 @@ contains | |
| @:ACC_SETUP_SFs(corrected_gps) | ||
|
|
||
| $:GPU_ENTER_DATA(copyin='[num_gps]') | ||
| $:GPU_UPDATE(device='[ib_second_order_vel, ib_ip_min_dist]') | ||
|
|
||
| if (collision_model > 0) call s_initialize_collisions_module() | ||
|
|
||
|
|
@@ -195,7 +196,7 @@ contains | |
| real(wp), intent(in) :: rho, pres_IP | ||
| real(wp), intent(out) :: pres_GP | ||
|
|
||
| pres_GP = pres_IP/min(max(1._wp - (abs(gp%levelset) + gp%ip_dist) & | ||
| & *rho/pres_IP*dot_product(patch_ib(gp_patch_id)%force/patch_ib(gp_patch_id)%mass, & | ||
| & gp%levelset_norm), 5.e-1_wp), 2._wp) | ||
|
|
||
|
|
@@ -259,6 +260,9 @@ contains | |
| if (buf > 0._wp) vel_GP = vel_GP + v_blow_eff*norm/buf | ||
| end if | ||
|
|
||
| ! extrapolate linearly through the boundary intercept, Mittal et al. (2008) | ||
| if (ib_second_order_vel) vel_GP = vel_GP - abs(gp%levelset)/max(gp%ip_dist, sgm_eps)*(vel_IP - vel_GP) | ||
|
There was a problem hiding this comment.
When AGENTS.md reference: AGENTS.md:L5-L5 Useful? React with 👍 / 👎. |
||
|
|
||
| end subroutine s_compute_ghost_point_velocity | ||
|
|
||
| !> Update the conservative variables at the ghost points | ||
|
|
@@ -708,6 +712,7 @@ contains | |
| impure subroutine s_compute_image_points() | ||
|
|
||
| real(wp) :: dist | ||
| real(wp) :: min_cell_width | ||
| real(wp), dimension(3) :: norm | ||
| real(wp), dimension(3) :: physical_loc | ||
| real(wp) :: temp_loc | ||
|
|
@@ -723,8 +728,8 @@ contains | |
|
|
||
| bounds_error = .false. | ||
|
|
||
| $:GPU_PARALLEL_LOOP(private='[q, gp, i, j, k, physical_loc, patch_id, dist, min_cell_width, norm, dim, bound, dir, index, & | ||
| & temp_loc, s_cc]', copy='[bounds_error]', present='[ghost_points]') | ||
| do q = 1, num_gps | ||
| gp = ghost_points(q) | ||
| i = gp%loc(1) | ||
|
|
@@ -742,7 +747,14 @@ contains | |
| patch_id = gp%ib_patch_id | ||
| dist = abs(real(gp%levelset, kind=wp)) | ||
| norm(:) = gp%levelset_norm | ||
| ghost_points(q)%ip_dist = dist | ||
| ! keep the image point stencil off the wall | ||
| if (ib_second_order_vel) then | ||
| min_cell_width = min(dx(i), dy(j)) | ||
| if (p > 0) min_cell_width = min(min_cell_width, dz(k)) | ||
| ghost_points(q)%ip_dist = max(dist, ib_ip_min_dist*min_cell_width) | ||
|
There was a problem hiding this comment.
When AGENTS.md reference: AGENTS.md:L5-L5 Useful? React with 👍 / 👎. |
||
| end if | ||
| ghost_points(q)%ip_loc(:) = physical_loc(:) + (dist + ghost_points(q)%ip_dist)*norm(:) | ||
|
|
||
| ! Find the closest grid point to the image point | ||
| do dim = 1, num_dims | ||
|
|
@@ -933,16 +945,12 @@ contains | |
| subroutine s_compute_interpolation_coeffs() | ||
|
|
||
| real(wp), dimension(2, 2, 2) :: dist | ||
| real(wp), dimension(2, 2, 2) :: fluid_mask | ||
| real(wp), dimension(2, 2, 2) :: interp_coeffs | ||
| type(ghost_point) :: gp | ||
| integer :: q, i, j, k, ii, jj, kk !< Grid indexes and iterators | ||
|
|
||
| $:GPU_PARALLEL_LOOP(private='[q, i, j, k, ii, jj, kk, dist, gp, interp_coeffs, fluid_mask]', present='[ghost_points]') | ||
| do q = 1, num_gps | ||
| gp = ghost_points(q) | ||
| ! Get the interpolation points | ||
|
|
@@ -956,7 +964,6 @@ contains | |
|
|
||
| ! get the distance to a cell in each direction | ||
| dist = 0._wp | ||
| do ii = 0, 1 | ||
| do jj = 0, 1 | ||
| if (p == 0) then | ||
|
|
@@ -971,51 +978,17 @@ contains | |
| end do | ||
| end do | ||
|
|
||
| ! flag the fluid cells of the stencil | ||
| fluid_mask = 0._wp | ||
| do kk = 0, min(p, 1) | ||
| do jj = 0, 1 | ||
| do ii = 0, 1 | ||
| if (ib_markers%sf(i + ii, j + jj, k + kk) == 0) fluid_mask(ii + 1, jj + 1, kk + 1) = 1._wp | ||
| end do | ||
| end do | ||
| end do | ||
|
|
||
| call s_interpolation_weights(gp, i, j, k, dist, fluid_mask, interp_coeffs) | ||
|
|
||
| ! An image point buried in a neighboring IB gets no weights at all: its ghost point is | ||
| ! averaged from its own neighbors instead, in s_ibm_correct_state. | ||
|
|
@@ -1027,6 +1000,53 @@ contains | |
|
|
||
| end subroutine s_compute_interpolation_coeffs | ||
|
|
||
| !> Fluid-masked image-point weights: bilinear for second-order IBM, else inverse distance squared | ||
| subroutine s_interpolation_weights(gp, i, j, k, dist, fluid_mask, interp_coeffs) | ||
|
|
||
| $:GPU_ROUTINE(parallelism='[seq]') | ||
|
|
||
| type(ghost_point), intent(in) :: gp | ||
| integer, intent(in) :: i, j, k !< lower corner of the stencil | ||
| real(wp), dimension(2, 2, 2), intent(in) :: dist, fluid_mask | ||
| real(wp), dimension(2, 2, 2), intent(out) :: interp_coeffs | ||
| real(wp), dimension(2, 2, 2) :: weights | ||
| real(wp), dimension(2, 3) :: axis_weight !< (stencil side, dimension) linear weights | ||
| real(wp) :: weight_sum | ||
| integer :: ii, jj, kk, num_z_layers | ||
|
|
||
| num_z_layers = min(p, 1) + 1 | ||
|
|
||
| if (ib_second_order_vel) then | ||
| ! bilinear weights, Mittal et al. (2008) | ||
| axis_weight(2, 1) = (gp%ip_loc(1) - x_cc(i))/(x_cc(i + 1) - x_cc(i)) | ||
| axis_weight(2, 2) = (gp%ip_loc(2) - y_cc(j))/(y_cc(j + 1) - y_cc(j)) | ||
| axis_weight(2, 3) = 0._wp | ||
| if (p /= 0) axis_weight(2, 3) = (gp%ip_loc(3) - z_cc(k))/(z_cc(k + 1) - z_cc(k)) | ||
| axis_weight(2,:) = min(max(axis_weight(2,:), 0._wp), 1._wp) ! clamp to the cell edges | ||
| axis_weight(1,:) = 1._wp - axis_weight(2,:) | ||
| do kk = 1, num_z_layers | ||
| do jj = 1, 2 | ||
| do ii = 1, 2 | ||
| weights(ii, jj, kk) = axis_weight(ii, 1)*axis_weight(jj, 2)*axis_weight(kk, 3) | ||
| end do | ||
| end do | ||
| end do | ||
| else | ||
| ! inverse distance squared, Chaudhuri et al. JCP (2011) | ||
| if (minval(dist(:,:,1:num_z_layers)) <= 1.e-16_wp) then | ||
| weights(:,:,1:num_z_layers) = merge(1._wp, 0._wp, dist(:,:,1:num_z_layers) <= 1.e-16_wp) | ||
| else | ||
| weights(:,:,1:num_z_layers) = 1._wp/dist(:,:,1:num_z_layers)**2 | ||
| end if | ||
| end if | ||
|
|
||
| interp_coeffs = 0._wp | ||
| weights(:,:,1:num_z_layers) = fluid_mask(:,:,1:num_z_layers)*weights(:,:,1:num_z_layers) | ||
| weight_sum = sum(weights(:,:,1:num_z_layers)) | ||
| if (weight_sum > 0._wp) interp_coeffs(:,:,1:num_z_layers) = weights(:,:,1:num_z_layers)/weight_sum | ||
|
|
||
| end subroutine s_interpolation_weights | ||
|
|
||
| !> Interpolate primitive variables to a ghost point's image point using bilinear or trilinear interpolation | ||
| subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, mv_IP, & | ||
| & nmom_IP, pb_in, mv_in, presb_IP, massv_IP, Ys_IP) | ||
|
|
||
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
For an image point that is not symmetrically placed among its stencil centers,
vel_IPcomes froms_compute_interpolation_coeffs, whose normalized1/dist**2weights at lines 1017-1028 do not reproduce even a linear velocity field. This new no-slip extrapolation therefore propagates an O(h) image-velocity error into the ghost value and cannot provide the advertised second-order surface velocity; use bilinear/trilinear linear-exact weights or an equivalent correction. The addedCircle -> second order velgoldens use only onen=49grid, so they cannot verify the claimed convergence order.AGENTS.md reference: AGENTS.md:L5-L5
Useful? React with 👍 / 👎.