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 diff --git a/src/common/m_compute_levelset.fpp b/src/common/m_compute_levelset.fpp index c7a3df21cf..da2b29956c 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_derived_types.fpp b/src/common/m_derived_types.fpp index f3ef91e831..5c36a7ae11 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 1bb00c44a3..3db86ceab5 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(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 @@ -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(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/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index ea679ba6f5..c95e1a2946 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -572,6 +572,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 0b8ec1ab81..7122102f50 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -95,6 +95,7 @@ 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 @@ -102,6 +103,7 @@ contains 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]') @@ -1150,6 +1152,60 @@ 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, num_cells_local + 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. & + patch_ib(ib_marker)%geometry == 5 .or. & + patch_ib(ib_marker)%geometry == 11 .or. & + patch_ib(ib_marker)%geometry == 12) then + + 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 + 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 + 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 + end do + + ! 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_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) + + ! 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 + 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_centroid_offset + subroutine s_compute_moment_of_inertia(ib_marker, axis) real(wp), dimension(3), intent(in) :: axis !< the axis about which we compute the moment. Only required in 3D. diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index b196b3afa3..6aaedae602 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: @@ -625,7 +626,10 @@ 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 + ! 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 call s_compute_ib_forces(q_prim_vf, 1._wp/fluid_pp(1)%Re(1)) diff --git a/toolchain/mfc/case.py b/toolchain/mfc/case.py index 269c8438bf..3125490dd5 100644 --- a/toolchain/mfc/case.py +++ b/toolchain/mfc/case.py @@ -7,12 +7,15 @@ from .state import ARG from .run import case_dicts - QPVF_IDX_VARS = { 'alpha_rho': 'contxb', 'vel' : 'momxb', 'pres': 'E_idx', '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)', '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 @dataclasses.dataclass(init=False) @@ -48,7 +51,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 @@ -95,8 +98,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]+\)%{re.escape(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 +196,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 @@ -266,9 +347,16 @@ 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) + return out + def get_fpp(self, target, print = True) -> str: def _prepend() -> str: return f"""\ diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 9240201527..dfdc13b38d 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -121,9 +121,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 @@ -372,9 +370,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