From 89a0ac5517aea9498685250dc5c4705bcac6840e Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Fri, 12 Dec 2025 13:48:00 -0500 Subject: [PATCH 01/11] Added ellipse marker patch --- src/common/m_ib_patches.fpp | 44 +++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index c0d95e6f12..18fcf5968e 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -119,6 +119,9 @@ contains ! STL+IBM patch elseif (patch_ib(i)%geometry == 5) then call s_ib_model(i, ib_markers_sf, levelset, levelset_norm) + elseif (patch_ib(i)%geometry == 6) then + call s_ib_ellipse(i, ib_markers_sf) + call s_ellipse_levelset(i, levelset, levelset_norm) end if end do !> @} @@ -762,6 +765,47 @@ contains end subroutine s_ib_cylinder + + subroutine s_ib_ellipse(patch_id, ib_markers_sf) + + integer, intent(in) :: patch_id + integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf + + integer :: i, j, k !< generic loop iterators + real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame + real(wp), dimension(1:2) :: length, center !< x and y coordinates in local IB frame + real(wp), dimension(1:3, 1:3) :: inverse_rotation + + ! Transferring the rectangle's centroid and length information + center(1) = patch_ib(patch_id)%x_centroid + center(2) = patch_ib(patch_id)%y_centroid + length(1) = patch_ib(patch_id)%length_x + length(2) = patch_ib(patch_id)%length_y + inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :) + + ! Checking whether the rectangle covers a particular cell in the + ! domain and verifying whether the current patch has the permission + ! to write to that cell. If both queries check out, the primitive + ! variables of the current patch are assigned to this cell. + $:GPU_PARALLEL_LOOP(private='[i,j, xy_local]', copy='[ib_markers_sf]',& + & copyin='[patch_id,center,length,inverse_rotation,x_cc,y_cc]', collapse=2) + do j = 0, n + do i = 0, m + ! get the x and y coordinates in the local IB frame + xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] + xy_local = matmul(inverse_rotation, xy_local) + + ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1 + if ((xy_local(1)/length(1))**2 + (xy_local(1)/length(1))**2) then + ! Updating the patch identities bookkeeping variable + ib_markers_sf(i, j, 0) = patch_id + end if + end do + end do + $:END_GPU_PARALLEL_LOOP() + + end subroutine s_ib_ellipse + !> The STL patch is a 2/3D geometry that is imported from an STL file. !! @param patch_id is the patch identifier !! @param ib_markers_sf Array to track patch ids From ef446f5addc81da2c82c13ec236ab4da5464d222 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Fri, 12 Dec 2025 16:40:40 -0500 Subject: [PATCH 02/11] Added ellipse IB patch --- src/common/m_compute_levelset.fpp | 59 +++++++++++++++++++++++++- src/common/m_ib_patches.fpp | 11 ++--- src/pre_process/m_check_ib_patches.fpp | 2 + 3 files changed, 66 insertions(+), 6 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index f343c0c47c..5346eb1b8b 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -23,7 +23,8 @@ module m_compute_levelset s_3D_airfoil_levelset, & s_rectangle_levelset, & s_cuboid_levelset, & - s_sphere_levelset + s_sphere_levelset, & + s_ellipse_levelset contains @@ -336,6 +337,62 @@ contains end subroutine s_rectangle_levelset + subroutine s_ellipse_levelset(ib_patch_id, levelset, levelset_norm) + + type(levelset_field), intent(INOUT), optional :: levelset + type(levelset_norm_field), intent(INOUT), optional :: levelset_norm + + integer, intent(in) :: ib_patch_id + real(wp) :: ellipse_coeffs(2) ! a and b in the ellipse equation + real(wp) :: quadratic_coeffs(3) ! A, B, C in the quadratic equation to compute levelset + + real(wp) :: length_x, length_y + real(wp), dimension(1:3) :: xy_local, normal_vector !< x and y coordinates in local IB frame + real(wp), dimension(2) :: center !< x and y coordinates in local IB frame + real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation + + integer :: i, j, k !< Loop index variables + integer :: idx !< Shortest path direction indicator + + length_x = patch_ib(ib_patch_id)%length_x + length_y = patch_ib(ib_patch_id)%length_y + center(1) = patch_ib(ib_patch_id)%x_centroid + center(2) = patch_ib(ib_patch_id)%y_centroid + inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) + rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) + + ellipse_coeffs(1) = 0.5_wp * length_x + ellipse_coeffs(2) = 0.5_wp * length_y + + $:GPU_PARALLEL_LOOP(private='[i,j,k,idx,quadratic_coeffs,xy_local,normal_vector]', & + & copyin='[ib_patch_id,center,ellipse_coeffs,inverse_rotation,rotation]', collapse=2) + do i = 0, m + do j = 0, n + xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] + xy_local = matmul(inverse_rotation, xy_local) + + ! we will get NaNs in the levelset if we compute this outside the ellipse + if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then + + normal_vector = xy_local + normal_vector(2) = normal_vector(2) * (ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transofmration method + normal_vector = normal_vector / sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector + levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame + + ! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3 + quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2 + quadratic_coeffs(2) = 2._wp * ((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2))) + quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp + + ! compute the levelset with the quadratic equation + levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp * (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / quadratic_coeffs(1) + end if + end do + end do + $:END_GPU_PARALLEL_LOOP() + + end subroutine s_ellipse_levelset + subroutine s_cuboid_levelset(ib_patch_id, levelset, levelset_norm) type(levelset_field), intent(INOUT), optional :: levelset diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index 18fcf5968e..8468d0b859 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -773,14 +773,15 @@ contains integer :: i, j, k !< generic loop iterators real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame - real(wp), dimension(1:2) :: length, center !< x and y coordinates in local IB frame + real(wp), dimension(1:2) :: ellipse_coeffs !< a and b in the ellipse coefficients + real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame real(wp), dimension(1:3, 1:3) :: inverse_rotation ! Transferring the rectangle's centroid and length information center(1) = patch_ib(patch_id)%x_centroid center(2) = patch_ib(patch_id)%y_centroid - length(1) = patch_ib(patch_id)%length_x - length(2) = patch_ib(patch_id)%length_y + ellipse_coeffs(1) = 0.5_wp * patch_ib(patch_id)%length_x + ellipse_coeffs(2) = 0.5_wp * patch_ib(patch_id)%length_y inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :) ! Checking whether the rectangle covers a particular cell in the @@ -788,7 +789,7 @@ contains ! to write to that cell. If both queries check out, the primitive ! variables of the current patch are assigned to this cell. $:GPU_PARALLEL_LOOP(private='[i,j, xy_local]', copy='[ib_markers_sf]',& - & copyin='[patch_id,center,length,inverse_rotation,x_cc,y_cc]', collapse=2) + & copyin='[patch_id,center,ellipse_coeffs,inverse_rotation,x_cc,y_cc]', collapse=2) do j = 0, n do i = 0, m ! get the x and y coordinates in the local IB frame @@ -796,7 +797,7 @@ contains xy_local = matmul(inverse_rotation, xy_local) ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1 - if ((xy_local(1)/length(1))**2 + (xy_local(1)/length(1))**2) then + if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then ! Updating the patch identities bookkeeping variable ib_markers_sf(i, j, 0) = patch_id end if diff --git a/src/pre_process/m_check_ib_patches.fpp b/src/pre_process/m_check_ib_patches.fpp index 7e1f9aa1d7..19493814e5 100644 --- a/src/pre_process/m_check_ib_patches.fpp +++ b/src/pre_process/m_check_ib_patches.fpp @@ -62,6 +62,8 @@ contains else if (patch_ib(i)%geometry == 5 .or. & patch_ib(i)%geometry == 12) then call s_check_model_ib_patch_geometry(i) + else if (patch_ib(i)%geometry == 6) then + print *, "Ellipse Patch" else call s_prohibit_abort("Invalid IB patch", & "patch_ib("//trim(iStr)//")%geometry must be "// & From 2f9b1beac4c083033fc8475f4568725683ed573f Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Fri, 12 Dec 2025 16:51:26 -0500 Subject: [PATCH 03/11] Added an ellipse immersed boundary patch --- examples/2D_ibm_ellipse/case.py | 97 +++++++++++++++++++++++++++++++ src/common/m_compute_levelset.fpp | 2 +- 2 files changed, 98 insertions(+), 1 deletion(-) create mode 100644 examples/2D_ibm_ellipse/case.py diff --git a/examples/2D_ibm_ellipse/case.py b/examples/2D_ibm_ellipse/case.py new file mode 100644 index 0000000000..1c5d8052f7 --- /dev/null +++ b/examples/2D_ibm_ellipse/case.py @@ -0,0 +1,97 @@ +import json +import math + +Mu = 1.84e-05 +gam_a = 1.4 + +# Configuring case dictionary +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + # For these computations, the cylinder is placed at the (0,0,0) + # domain origin. + # axial direction + "x_domain%beg": 0.0e00, + "x_domain%end": 6.0e-03, + # r direction + "y_domain%beg": 0.0e00, + "y_domain%end": 6.0e-03, + "cyl_coord": "F", + "m": 250, + "n": 250, + "p": 0, + "dt": 0.5e-5, + "t_step_start": 0, + "t_step_stop": 1000, # 3000 + "t_step_save": 10, # 10 + # Simulation Algorithm Parameters + # Only one patches are necessary, the air tube + "num_patches": 1, + # Use the 5 equation model + "model_eqns": 2, + "alt_soundspeed": "F", + # One fluids: air + "num_fluids": 1, + # time step + "mpp_lim": "F", + # Correct errors when computing speed of sound + "mixture_err": "T", + # Use TVD RK3 for time marching + "time_stepper": 3, + # Use WENO5 + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_Re_flux": "T", + "weno_avg": "T", + "avg_state": 2, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "T", + "riemann_solver": 2, + "wave_speeds": 1, + # We use ghost-cell + "bc_x%beg": -3, + "bc_x%end": -3, + "bc_y%beg": -3, + "bc_y%end": -3, + # Set IB to True and add 1 patch + "ib": "T", + "num_ibs": 1, + "viscous": "T", + # Formatted Database Files Structure Parameters + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "E_wrt": "T", + "parallel_io": "T", + # Patch: Constant Tube filled with air + # Specify the cylindrical air tube grid geometry + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 3.0e-03, + # Uniform medium density, centroid is at the center of the domain + "patch_icpp(1)%y_centroid": 3.0e-03, + "patch_icpp(1)%length_x": 6.0e-03, + "patch_icpp(1)%length_y": 6.0e-03, + # Specify the patch primitive variables + "patch_icpp(1)%vel(1)": 0.05e00, + "patch_icpp(1)%vel(2)": 0.0e00, + "patch_icpp(1)%pres": 1.0e00, + "patch_icpp(1)%alpha_rho(1)": 1.0e00, + "patch_icpp(1)%alpha(1)": 1.0e00, + # Patch: Cylinder Immersed Boundary + "patch_ib(1)%geometry": 6, + "patch_ib(1)%x_centroid": 1.5e-03, + "patch_ib(1)%y_centroid": 3.0e-03, + "patch_ib(1)%length_x": 0.4e-03, + "patch_ib(1)%length_y": 0.2e-03, + "patch_ib(1)%slip": "F", + # Fluids Physical Parameters + "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00), # 2.50(Not 1.40) + "fluid_pp(1)%pi_inf": 0, + "fluid_pp(1)%Re(1)": 2500000, + } + ) +) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index 5346eb1b8b..69ddda5fd5 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -384,7 +384,7 @@ contains quadratic_coeffs(2) = 2._wp * ((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2))) quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp - ! compute the levelset with the quadratic equation + ! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp * (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / quadratic_coeffs(1) end if end do From d54bc5bd0b05c67c6305997945d2857ba85a8551 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 14 Jan 2026 11:22:44 -0500 Subject: [PATCH 04/11] initial implementation of analytic moving IB --- src/simulation/m_time_steppers.fpp | 4 +- toolchain/mfc/case.py | 90 +++++++++++++++++++++++++++++- 2 files changed, 92 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index f17689d6bc..420132ad83 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -625,7 +625,9 @@ contains patch_ib(i)%vel = (rk_coef(s, 1)*patch_ib(i)%step_vel + rk_coef(s, 2)*patch_ib(i)%vel)/rk_coef(s, 4) patch_ib(i)%angular_vel = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel + rk_coef(s, 2)*patch_ib(i)%angular_vel)/rk_coef(s, 4) - if (patch_ib(i)%moving_ibm == 2) then ! if we are using two-way coupling, apply force and torque + if (patch_ib(i)%moving_ibm == 1) then + @:mib_analytical() + else if (patch_ib(i)%moving_ibm == 2) then ! if we are using two-way coupling, apply force and torque ! compute the force and torque on the IB from the fluid call s_compute_ib_forces(q_prim_vf(E_idx)) diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 269c8438bf..5f7876ab0d 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -13,6 +13,10 @@ 'alpha': 'advxb', 'tau_e': 'stress_idx%beg', 'Y': 'chemxb', 'cf_val': 'c_idx', 'Bx': 'B_idx%beg', 'By': 'B_idx%end-1', 'Bz': 'B_idx%end', } + +MIBM_ANALYTIC_VARS = [ + 'vel(1)', 'vel(2)', 'vel(3)', 'angluar_vel(1)', 'angluar_vel(2)', 'angluar_vel(3)' +] # "B_idx%end - 1" not "B_idx%beg + 1" must be used because 1D does not have Bx @dataclasses.dataclass(init=False) @@ -95,8 +99,20 @@ def __is_ic_analytical(self, key: str, val: str) -> bool: return False + def __is_mib_analytical(self, key: str, val: str) -> bool: + '''Is this initial condition analytical? + More precisely, is this an arbitrary expression or a string representing a number?''' + if common.is_number(val) or not isinstance(val, str): + return False + + for variable in MIBM_ANALYTIC_VARS: + if re.match(fr'^patch_ib\([0-9]+\)%{variable}', key): + return True + + return False + # pylint: disable=too-many-locals - def __get_pre_fpp(self, print: bool) -> str: + def __get_analytic_ic_fpp(self, print: bool) -> str: # generates the content of an FFP file that will hold the functions for # some initial condition DATA = { @@ -181,6 +197,72 @@ def rhs_replace(match): #:def analytical() {f'{chr(10)}{chr(10)}'.join(srcs)} #:enddef +""" + return content + + # gets the analytic description of a moving IB's velocity and rotation rate + def __get_analytic_mib_fpp(self, print: bool) -> str: + # iterates over the parameters and checks if they are defined as an + # analytical function. If so, append it to the `patches`` object + ib_patches = {} + + for key, val in self.params.items(): + if not self.__is_mib_analytical(key, val): + continue + + patch_id = re.search(r'[0-9]+', key).group(0) + + if patch_id not in ib_patches: + ib_patches[patch_id] = [] + + ib_patches[patch_id].append((key, val)) + + srcs = [] + + # for each analytical patch that is required to be added, generate + # the string that contains that function. + for pid, items in ib_patches.items(): + + # function that defines how we will replace variable names with + # values from the case file + def rhs_replace(match): + return { + 'x': 'x_cc(i)', 'y': 'y_cc(j)', 'z': 'z_cc(k)', + 't': 'mytime', + + 'r': f'patch_ib({pid})%radius', + + 'e' : f'{math.e}', + }.get(match.group(), match.group()) + + lines = [] + # perform the replacement of strings for each analytic function + # to generate some fortran string representing the code passed in + for attribute, expr in items: + if print: + cons.print(f"* Codegen: {attribute} = {expr}") + + lhs = attribute + rhs = re.sub(r"[a-zA-Z]+", rhs_replace, expr) + + lines.append(f" {lhs} = {rhs}") + + # concatenates all of the analytic lines into a single string with + # each element separated by new line characters. Then write those + # new lines as a fully concatenated string with fortran syntax + srcs.append(f"""\ + if (i == {pid}) then +{f'{chr(10)}'.join(lines)} + end if\ +""") + + content = f"""\ +! This file was generated by MFC. It is only used when we analytically +! parameterize the velocity and rotation rate of a moving IB. + +#:def mib_analytical() +{f'{chr(10)}{chr(10)}'.join(srcs)} +#:enddef """ return content @@ -269,6 +351,12 @@ def __get_sim_fpp(self, print: bool) -> str: # We need to also include the pre_processing includes so that common subroutines have access to the @:analytical function return out + f"\n{self.__get_pre_fpp(print)}" + + def __get_pre_fpp(self, print: bool) -> str: + out = self.__get_analytic_ic_fpp(print) + out += self.__get_analytical_mib_fpp(print) + return out + def get_fpp(self, target, print = True) -> str: def _prepend() -> str: return f"""\ From 697b5a05e76f37ca8e5319f6699dc18972353dfd Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Wed, 14 Jan 2026 13:48:22 -0500 Subject: [PATCH 05/11] Works in simple example --- src/simulation/m_time_steppers.fpp | 2 ++ toolchain/mfc/case.py | 9 ++++++--- toolchain/mfc/run/case_dicts.py | 6 ++---- 3 files changed, 10 insertions(+), 7 deletions(-) diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 420132ad83..10cd51dcd4 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -3,6 +3,7 @@ !! @brief Contains module m_time_steppers #:include 'macros.fpp' +#:include 'case.fpp' !> @brief The following module features a variety of time-stepping schemes. !! Currently, it includes the following Runge-Kutta (RK) algorithms: @@ -626,6 +627,7 @@ contains patch_ib(i)%angular_vel = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel + rk_coef(s, 2)*patch_ib(i)%angular_vel)/rk_coef(s, 4) if (patch_ib(i)%moving_ibm == 1) then + ! plug in analytic velocities for 1-way coupling, if it exists @:mib_analytical() else if (patch_ib(i)%moving_ibm == 2) then ! if we are using two-way coupling, apply force and torque ! compute the force and torque on the IB from the fluid diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 5f7876ab0d..28e6fa0931 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -7,6 +7,8 @@ from .state import ARG from .run import case_dicts +from pprint import pprint + QPVF_IDX_VARS = { 'alpha_rho': 'contxb', 'vel' : 'momxb', 'pres': 'E_idx', @@ -52,7 +54,7 @@ def get_inp(self, _target) -> str: dict_str = "" for key, val in self.params.items(): if key in MASTER_KEYS and key not in case_dicts.IGNORE: - if self.__is_ic_analytical(key, val): + if self.__is_ic_analytical(key, val) or self.__is_mib_analytical(key, val): dict_str += f"{key} = 0d0\n" ignored.append(key) continue @@ -106,7 +108,7 @@ def __is_mib_analytical(self, key: str, val: str) -> bool: return False for variable in MIBM_ANALYTIC_VARS: - if re.match(fr'^patch_ib\([0-9]+\)%{variable}', key): + if re.match(fr'^patch_ib\([0-9]+\)%{re.escape(variable)}', key): return True return False @@ -348,13 +350,14 @@ def __get_sim_fpp(self, print: bool) -> str: else: out = "" + out = out + self.__get_analytic_mib_fpp(print) + # We need to also include the pre_processing includes so that common subroutines have access to the @:analytical function return out + f"\n{self.__get_pre_fpp(print)}" def __get_pre_fpp(self, print: bool) -> str: out = self.__get_analytic_ic_fpp(print) - out += self.__get_analytical_mib_fpp(print) return out def get_fpp(self, target, print = True) -> str: diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index a63af7267f..605f63ae61 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -119,9 +119,7 @@ def analytic(self): PRE_PROCESS[f"patch_ib({ib_id})%{real_attr}"] = ty for dir_id in range(1, 4): - PRE_PROCESS[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL PRE_PROCESS[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL - PRE_PROCESS[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL for cmp_id, cmp in enumerate(["x", "y", "z"]): cmp_id += 1 @@ -368,9 +366,9 @@ def analytic(self): SIMULATION[f"patch_ib({ib_id})%{real_attr}"] = ty for dir_id in range(1, 4): - SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL + SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic() SIMULATION[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL - SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL + SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL.analytic() for cmp_id, cmp in enumerate(["x", "y", "z"]): cmp_id += 1 From 6ed835071b1cd305d3c3028a0648547b67198042 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Thu, 15 Jan 2026 11:43:53 -0500 Subject: [PATCH 06/11] Added beginning of centroid offset for calculations --- src/common/m_derived_types.fpp | 1 + src/common/m_ib_patches.fpp | 2 + src/pre_process/m_global_parameters.fpp | 1 + src/simulation/m_ibm.fpp | 53 ++++++++++++++++++++++++- 4 files changed, 56 insertions(+), 1 deletion(-) diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 1d82405c64..da278a2b2b 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -296,6 +296,7 @@ module m_derived_types !! is specified through its x-, y- and z-coordinates, respectively. real(wp) :: step_x_centroid, step_y_centroid, step_z_centroid !< !! Centroid locations of intermediate steps in the time_stepper module + real(wp), dimension(1:3) :: centroid_offset ! offset of center of mass from computed cell center for odd-shaped IBs real(wp), dimension(1:3) :: angles real(wp), dimension(1:3) :: step_angles diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index 8468d0b859..a271e62323 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -279,6 +279,7 @@ contains do i = 0, m xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates + xy_local = xy_local + patch_ib(ib_marker)%centroid_offset ! airfoils are a patch that require a centroid offset if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then xa = xy_local(1)/ca_in @@ -433,6 +434,7 @@ contains do i = 0, m xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates + xyz_local = xyz_local + patch_ib(ib_marker)%centroid_offset ! airfoils are a patch that require a centroid offset if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index dc2ed9ac71..2896b98bcf 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -562,6 +562,7 @@ contains patch_ib(i)%angular_vel(:) = 0._wp patch_ib(i)%mass = dflt_real patch_ib(i)%moment = dflt_real + patch_ib(i)%centroid_offset(:) = 0._wp ! sets values of a rotation matrix which can be used when calculating rotations patch_ib(i)%rotation_matrix = 0._wp diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 048032b9ae..d972d12092 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -93,13 +93,14 @@ contains integer :: i, j, k integer :: max_num_gps, max_num_inner_gps + ! do all set up for moving immersed boundaries moving_immersed_boundary_flag = .false. do i = 1, num_ibs if (patch_ib(i)%moving_ibm /= 0) then moving_immersed_boundary_flag = .true. - end if call s_update_ib_rotation_matrix(i) + call s_compute_centroid_offset(i) end do $:GPU_ENTER_DATA(copyin='[patch_ib]') @@ -1080,6 +1081,56 @@ contains end subroutine s_finalize_ibm_module + !> Computes the center of mass for IB patch types where we are unable to determine their center of mass analytically. + !> These patches include things like NACA airfoils and STL models + subroutine s_compute_centroid_offset(ib_marker) + + integer, intent(in) :: ib_marker + + integer :: i, j, k, num_cells + real(wp), dimension(1:3) :: center_of_mass + + ! Offset only needs to be computes for specific geometries + if (patch_ib(ib_marker)%geometry == 4 .or. & + patch_ib(ib_marker)%geometry == 5 .or. & + patch_ib(ib_marker)%geometry == 11 .or. & + patch_ib(ib_marker)%geometry == 12 .or. ) then + + center_of_mass = [0._wp, 0._wp, 0._wp] + num_cells = 0 + + ! get the summed mass distribution and number of cells to divide by + do i = 0, m + do j = 0, n + do k = 0, p + if (ib_markers%sf(i, j, k) == ib_marker) then + center_of_mass = center_of_mass + [x_cc(i), y_cc(j), z_cc(k)] + num_cells = num_cells + 1 + end if + end do + end do + end do + + ! reduce the mass contribution over all MPI ranks and compute COM + call s_mpi_allreduce_vectors_sum(center_of_mass, center_of_mass, 1, 3) + call s_mpi_allreduce_integer_sum(num_cells, num_cells) + center_of_mass = center_of_mass / real(num_cells, wp) ! compute the actual center of mass + + ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the current centroid + patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%_centroid] & + - center_of_mass + patch_ib(ib_marker)%x_centroid = center_of_mass(1) + patch_ib(ib_marker)%y_centroid = center_of_mass(2) + patch_ib(ib_marker)%z_centroid = center_of_mass(3) + + ! rotate the centroid offset back into the local coords of the IB + patch_ib(ib_marker)%centroid_offset = matmul(patch_ib(ib_marker)%rotation_matrix_inverse, patch_ib(ib_marker)%centroid_offset) + else + patch_ib(ib_marker)%centroid_offset(:) = [0._wp, 0._wp, 0._wp] + end if + + end subroutine s_compute_center_of_mass_offset + subroutine s_compute_moment_of_inertia(ib_marker, axis) real(wp), dimension(3), optional :: axis !< the axis about which we compute the moment. Only required in 3D. From 7bd28cd32630d59b66391b085834a31ad34e7678 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Thu, 15 Jan 2026 14:52:43 -0500 Subject: [PATCH 07/11] Works with pitching airfoil --- src/common/m_compute_levelset.fpp | 2 ++ src/common/m_ib_patches.fpp | 4 ++-- src/simulation/m_ibm.fpp | 24 ++++++++++++++---------- toolchain/mfc/case.py | 5 +---- 4 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index 69ddda5fd5..47f64c08ad 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -94,6 +94,7 @@ contains do j = 0, n xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate + xy_local = xy_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset if (xy_local(2) >= 0._wp) then ! finds the location on the airfoil grid with the minimum distance (closest) @@ -189,6 +190,7 @@ contains xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates + xyz_local = xyz_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset if (xyz_local(2) >= center(2)) then do k = 1, Np diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index a271e62323..38d0d6bfed 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -279,7 +279,7 @@ contains do i = 0, m xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates - xy_local = xy_local + patch_ib(ib_marker)%centroid_offset ! airfoils are a patch that require a centroid offset + xy_local = xy_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then xa = xy_local(1)/ca_in @@ -434,7 +434,7 @@ contains do i = 0, m xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates - xyz_local = xyz_local + patch_ib(ib_marker)%centroid_offset ! airfoils are a patch that require a centroid offset + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index d972d12092..941a5baad4 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1087,37 +1087,41 @@ contains integer, intent(in) :: ib_marker - integer :: i, j, k, num_cells + integer :: i, j, k, num_cells, num_cells_local real(wp), dimension(1:3) :: center_of_mass ! Offset only needs to be computes for specific geometries if (patch_ib(ib_marker)%geometry == 4 .or. & patch_ib(ib_marker)%geometry == 5 .or. & patch_ib(ib_marker)%geometry == 11 .or. & - patch_ib(ib_marker)%geometry == 12 .or. ) then + patch_ib(ib_marker)%geometry == 12) then center_of_mass = [0._wp, 0._wp, 0._wp] - num_cells = 0 + num_cells_local = 0 ! get the summed mass distribution and number of cells to divide by do i = 0, m do j = 0, n do k = 0, p if (ib_markers%sf(i, j, k) == ib_marker) then - center_of_mass = center_of_mass + [x_cc(i), y_cc(j), z_cc(k)] - num_cells = num_cells + 1 + num_cells_local = num_cells_local + 1 + center_of_mass = center_of_mass + [x_cc(i), y_cc(j), 0._wp] + if (num_dims == 3) center_of_mass(3) = center_of_mass(3) + z_cc(k) end if end do end do end do ! reduce the mass contribution over all MPI ranks and compute COM - call s_mpi_allreduce_vectors_sum(center_of_mass, center_of_mass, 1, 3) - call s_mpi_allreduce_integer_sum(num_cells, num_cells) - center_of_mass = center_of_mass / real(num_cells, wp) ! compute the actual center of mass + ! print *, "Before reduction ", center_of_mass, num_cells_local + call s_mpi_allreduce_sum(center_of_mass(1), center_of_mass(1)) + call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2)) + call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3)) + call s_mpi_allreduce_integer_sum(num_cells_local, num_cells) + center_of_mass = center_of_mass / real(num_cells_local, wp) ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the current centroid - patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%_centroid] & + patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid] & - center_of_mass patch_ib(ib_marker)%x_centroid = center_of_mass(1) patch_ib(ib_marker)%y_centroid = center_of_mass(2) @@ -1129,7 +1133,7 @@ contains patch_ib(ib_marker)%centroid_offset(:) = [0._wp, 0._wp, 0._wp] end if - end subroutine s_compute_center_of_mass_offset + end subroutine s_compute_centroid_offset subroutine s_compute_moment_of_inertia(ib_marker, axis) diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 28e6fa0931..3125490dd5 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -7,9 +7,6 @@ from .state import ARG from .run import case_dicts -from pprint import pprint - - QPVF_IDX_VARS = { 'alpha_rho': 'contxb', 'vel' : 'momxb', 'pres': 'E_idx', 'alpha': 'advxb', 'tau_e': 'stress_idx%beg', 'Y': 'chemxb', @@ -17,7 +14,7 @@ } MIBM_ANALYTIC_VARS = [ - 'vel(1)', 'vel(2)', 'vel(3)', 'angluar_vel(1)', 'angluar_vel(2)', 'angluar_vel(3)' + 'vel(1)', 'vel(2)', 'vel(3)', 'angular_vel(1)', 'angular_vel(2)', 'angular_vel(3)' ] # "B_idx%end - 1" not "B_idx%beg + 1" must be used because 1D does not have Bx From f40422a18251b2d67d00316ec3b649f610181b00 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Thu, 15 Jan 2026 15:09:15 -0500 Subject: [PATCH 08/11] Spelling and formatting --- src/common/m_compute_levelset.fpp | 12 ++++++------ src/common/m_ib_patches.fpp | 5 ++--- src/simulation/m_ibm.fpp | 8 ++++---- 3 files changed, 12 insertions(+), 13 deletions(-) diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index 47f64c08ad..da2b29956c 100644 --- a/src/common/m_compute_levelset.fpp +++ b/src/common/m_compute_levelset.fpp @@ -363,8 +363,8 @@ contains inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :) rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :) - ellipse_coeffs(1) = 0.5_wp * length_x - ellipse_coeffs(2) = 0.5_wp * length_y + ellipse_coeffs(1) = 0.5_wp*length_x + ellipse_coeffs(2) = 0.5_wp*length_y $:GPU_PARALLEL_LOOP(private='[i,j,k,idx,quadratic_coeffs,xy_local,normal_vector]', & & copyin='[ib_patch_id,center,ellipse_coeffs,inverse_rotation,rotation]', collapse=2) @@ -377,17 +377,17 @@ contains if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then normal_vector = xy_local - normal_vector(2) = normal_vector(2) * (ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transofmration method - normal_vector = normal_vector / sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector + normal_vector(2) = normal_vector(2)*(ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transformation method + normal_vector = normal_vector/sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame ! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3 quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2 - quadratic_coeffs(2) = 2._wp * ((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2))) + quadratic_coeffs(2) = 2._wp*((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2))) quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp ! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A - levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp * (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / quadratic_coeffs(1) + levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp*(-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp*quadratic_coeffs(1)*quadratic_coeffs(3)))/quadratic_coeffs(1) end if end do end do diff --git a/src/common/m_ib_patches.fpp b/src/common/m_ib_patches.fpp index 38d0d6bfed..dff617ef68 100644 --- a/src/common/m_ib_patches.fpp +++ b/src/common/m_ib_patches.fpp @@ -767,7 +767,6 @@ contains end subroutine s_ib_cylinder - subroutine s_ib_ellipse(patch_id, ib_markers_sf) integer, intent(in) :: patch_id @@ -782,8 +781,8 @@ contains ! Transferring the rectangle's centroid and length information center(1) = patch_ib(patch_id)%x_centroid center(2) = patch_ib(patch_id)%y_centroid - ellipse_coeffs(1) = 0.5_wp * patch_ib(patch_id)%length_x - ellipse_coeffs(2) = 0.5_wp * patch_ib(patch_id)%length_y + ellipse_coeffs(1) = 0.5_wp*patch_ib(patch_id)%length_x + ellipse_coeffs(2) = 0.5_wp*patch_ib(patch_id)%length_y inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :) ! Checking whether the rectangle covers a particular cell in the diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 941a5baad4..f0b8ab7eb3 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1098,13 +1098,13 @@ contains center_of_mass = [0._wp, 0._wp, 0._wp] num_cells_local = 0 - + ! get the summed mass distribution and number of cells to divide by do i = 0, m do j = 0, n do k = 0, p if (ib_markers%sf(i, j, k) == ib_marker) then - num_cells_local = num_cells_local + 1 + num_cells_local = num_cells_local + 1 center_of_mass = center_of_mass + [x_cc(i), y_cc(j), 0._wp] if (num_dims == 3) center_of_mass(3) = center_of_mass(3) + z_cc(k) end if @@ -1118,11 +1118,11 @@ contains call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2)) call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3)) call s_mpi_allreduce_integer_sum(num_cells_local, num_cells) - center_of_mass = center_of_mass / real(num_cells_local, wp) + center_of_mass = center_of_mass/real(num_cells_local, wp) ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the current centroid patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid] & - - center_of_mass + - center_of_mass patch_ib(ib_marker)%x_centroid = center_of_mass(1) patch_ib(ib_marker)%y_centroid = center_of_mass(2) patch_ib(ib_marker)%z_centroid = center_of_mass(3) From 541e727d70e01e303ab3af1b9b820ba05447c492 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Thu, 15 Jan 2026 15:10:09 -0500 Subject: [PATCH 09/11] Hey, the AI finally caught a bug --- src/simulation/m_ibm.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index f0b8ab7eb3..7491d1768e 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1118,7 +1118,7 @@ contains call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2)) call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3)) call s_mpi_allreduce_integer_sum(num_cells_local, num_cells) - center_of_mass = center_of_mass/real(num_cells_local, wp) + center_of_mass = center_of_mass/real(num_cells, wp) ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the current centroid patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid] & From cb7b845055e86bf07e17e1c9bff830d5884fed6e Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Thu, 15 Jan 2026 16:12:17 -0500 Subject: [PATCH 10/11] Added empty analytic function that can be found in documentation --- src/common/include/case.fpp | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/common/include/case.fpp b/src/common/include/case.fpp index 84d1d2cc14..e9d05e9edf 100644 --- a/src/common/include/case.fpp +++ b/src/common/include/case.fpp @@ -6,3 +6,8 @@ #:def analytical() #:enddef + +! For moving immersed boundaries in simulation +#:def mib_analytical() + +#:enddef \ No newline at end of file From 37f27304fbc236747fbd1c9f8bc5a5e0cf28afd6 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Fri, 16 Jan 2026 17:59:20 -0500 Subject: [PATCH 11/11] Fixed multi-rank issue on frontier --- src/simulation/m_ibm.fpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 737af54c14..7122102f50 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1159,7 +1159,7 @@ contains integer, intent(in) :: ib_marker integer :: i, j, k, num_cells, num_cells_local - real(wp), dimension(1:3) :: center_of_mass + real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local ! Offset only needs to be computes for specific geometries if (patch_ib(ib_marker)%geometry == 4 .or. & @@ -1167,7 +1167,7 @@ contains patch_ib(ib_marker)%geometry == 11 .or. & patch_ib(ib_marker)%geometry == 12) then - center_of_mass = [0._wp, 0._wp, 0._wp] + center_of_mass_local = [0._wp, 0._wp, 0._wp] num_cells_local = 0 ! get the summed mass distribution and number of cells to divide by @@ -1176,8 +1176,8 @@ contains do k = 0, p if (ib_markers%sf(i, j, k) == ib_marker) then num_cells_local = num_cells_local + 1 - center_of_mass = center_of_mass + [x_cc(i), y_cc(j), 0._wp] - if (num_dims == 3) center_of_mass(3) = center_of_mass(3) + z_cc(k) + center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp] + if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k) end if end do end do @@ -1185,9 +1185,9 @@ contains ! reduce the mass contribution over all MPI ranks and compute COM ! print *, "Before reduction ", center_of_mass, num_cells_local - call s_mpi_allreduce_sum(center_of_mass(1), center_of_mass(1)) - call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2)) - call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3)) + call s_mpi_allreduce_sum(center_of_mass_local(1), center_of_mass(1)) + call s_mpi_allreduce_sum(center_of_mass_local(2), center_of_mass(2)) + call s_mpi_allreduce_sum(center_of_mass_local(3), center_of_mass(3)) call s_mpi_allreduce_integer_sum(num_cells_local, num_cells) center_of_mass = center_of_mass/real(num_cells, wp)