diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 259b0bed4..c14158c67 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -773,6 +773,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca | `heat_ratio_wrt` | Logical | Add the specific heat ratio to the database | | `ib_force_wrt` | Logical | Record the immersed-boundary force history to `D/ib_forces.dat` (default off) | | `ib_force_stride` | Integer | Stride, in time steps, of the per-step immersed-boundary force record (default 1) | +| `ib_surface_wrt` | Logical | Write the wall temperature and gasified mass flux at each thermal/reacting IB surface point at every save (default off) | | `ib_state_wrt` | Logical | Parameter to handle writing IB state on saves and outputting the state as a point mesh to SILO files. | | `pi_inf_wrt` | Logical | Add the liquid stiffness function to the database | | `pres_inf_wrt` | Logical | Add the liquid stiffness to the formatted database | @@ -844,6 +845,8 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu - `ib_force_wrt` records the force, torque and kinematics of every immersed boundary in a single shared text file, `D/ib_forces.dat`, described below. It is off by default: the history is written every step, which at large rank counts is a cost a run should opt into rather than inherit. `ib_force_stride` writes only every N-th step, for runs long enough that the history itself becomes large. +- `ib_surface_wrt` writes, at every save, one text file per rank, `D/ib_surface__.dat`, with a line per surface point of each chemistry IB (`thermal_bc` /= 0 or `surface_reaction` = 1) and no header, columns `x y z ib nx ny nz area T_wall mdot`. `ib` is the global IB index, `(nx, ny, nz)` the level-set normal, `T_wall` the surface temperature [K] and `mdot` the gasified mass flux [kg/m²/s] from the surface solve (0 on an inert surface). `area` [m², or m per unit depth in 2D] is the surface each point stands for: summing `area`·`mdot` over an IB's lines gives its mass loss rate [kg/s], and `mdot`/ρ_solid is the local surface regression rate. The points are the ghost points within two cell sizes h = (cell volume)^(1/d) of the surface, each weighted by cell volume/(2h); for circles, spheres and cylinder sides the weight is also scaled by (R/(R − depth))^(d−1), the area ratio between the surface and the layer the point sits in. The areas then sum to the surface area without bias (a sphere at R = 12h: within 0.3%; a circle at R = 20h: within about 1%); other shapes keep the uncorrected band, which underestimates a convex surface by about (d−1)·h/R. + #### Immersed-boundary force history {#sec-ib-force-history} `D/ib_forces.dat` holds one fixed-width record per body per written step. Its twenty columns are diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 6157db87e..e05701731 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -518,6 +518,7 @@ contains ib_coefficient_of_friction = dflt_real ib_state_wrt = .false. ib_force_wrt = .false. + ib_surface_wrt = .false. ib_force_stride = 1 many_ib_patch_parallelism = .false. diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 614d0dc4b..905353df6 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -31,7 +31,8 @@ module m_ibm 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 - ; public :: s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module, s_report_ibm_surface + ; public :: s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module, s_report_ibm_surface, & + & s_write_ib_surface !> Ghost points at which the reacting-surface Newton solve did not reach its tolerance, so the point fell back to a chemically !! inert wall (keeping a prescribed Twall). Counted because that fallback is otherwise indistinguishable from a surface @@ -52,6 +53,13 @@ module m_ibm type(ghost_point), dimension(:), allocatable :: ghost_points $:GPU_DECLARE(create='[ghost_points]') + !> Surface record per ghost point from the latest s_ibm_correct_state, kept when ib_surface_wrt: (1) area weight, (2) wall + !! temperature, (3) gasified mass flux. Only ghost points within two cells of the surface carry a weight, so summing weight*flux + !! over them integrates over the surface; see s_record_gp_surface. + integer, parameter :: ib_surf_nvars = 3 + real(wp), allocatable, dimension(:,:) :: gp_surf + $:GPU_DECLARE(create='[gp_surf]') + integer :: num_gps !< Number of ghost points #if defined(MFC_OpenACC) $:GPU_DECLARE(create='[gp_layers, num_gps]') @@ -170,6 +178,11 @@ contains @:ALLOCATE(ghost_points(1:max_num_gps)) $:GPU_ENTER_DATA(copyin='[ghost_points]') + if (ib_surface_wrt) then + @:ALLOCATE(gp_surf(ib_surf_nvars, 1:max_num_gps)) + gp_surf = 0._wp + $:GPU_UPDATE(device='[gp_surf]') + end if ! Ghost-cell IBM, Tseng & Ferziger JCP (2003), Mittal & Iaccarino ARFM (2005) call s_find_ghost_points() call s_apply_levelset(ghost_points, num_gps) @@ -347,6 +360,7 @@ contains & surface_converged, vel_sum_g, E_ghost, alpha_q, alpha_rho_q, e_q]', & & reduction='[[n_not_converged, n_ill_posed]]', reductionOp='[+]', present='[ghost_points]') do i = 1, num_gps + if (ib_surface_wrt) gp_surf(:,i) = 0._wp gp = ghost_points(i) if (.not. gp%interp_valid) cycle j = gp%loc(1) @@ -434,6 +448,8 @@ contains if (patch_ib(patch_id)%thermal_bc == 1) T_s = patch_ib(patch_id)%Twall end if + if (ib_surface_wrt) call s_record_gp_surface(i, gp, T_s, mdot_s, surface_converged) + call s_blend_ghost_state(T_IP, T_s, Ys_IP, Ys_s, T_g, Ys_g) call get_mixture_molecular_weight(Ys_g, mw_g) @@ -702,6 +718,84 @@ contains end subroutine s_report_ibm_surface + !> Store ghost point i's surface record in gp_surf. The weight is the surface area the point stands for. Ghost points within two + !! cell sizes h = dV^(1/d) of the surface fill a band of volume ~2hA, so dV/(2h) each sums to the area A (a length in 2D) + !! whatever the surface's orientation to the grid. The band lies inside the solid, where a layer at depth s has area A (1 - + !! s/R)^(d-1) on a curved surface; for circles, spheres and cylinders the weight is scaled back by (R/(R - s))^(d-1), which + !! removes that bias (sampled on random lattice offsets: < 0.1% mean, 0.3% spread for a sphere at R = 12h). + subroutine s_record_gp_surface(i, gp, T_s, mdot_s, reacting) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: i + type(ghost_point), intent(in) :: gp + real(wp), intent(in) :: T_s, mdot_s + logical, intent(in) :: reacting + real(wp) :: d, dV, h, R + + d = abs(real(gp%levelset, kind=wp)) + dV = dx(gp%loc(1))*dy(gp%loc(2)) + if (num_dims == 3) dV = dV*dz(gp%loc(3)) + h = dV**(1._wp/real(num_dims, wp)) + if (.not. (d > 0._wp .and. d <= 2._wp*h)) return + + gp_surf(1, i) = dV/(2._wp*h) + R = patch_ib(gp%ib_patch_id)%radius + select case (patch_ib(gp%ib_patch_id)%geometry) + case (2, 8) + gp_surf(1, i) = gp_surf(1, i)*(R/(R - d))**(num_dims - 1) + case (10) ! curved side only; the flat caps need no correction + if (abs(gp%levelset_norm(f_cylinder_axis(patch_ib(gp%ib_patch_id)))) < 0.5_wp) gp_surf(1, i) = gp_surf(1, i)*R/(R - d) + end select + + gp_surf(2, i) = T_s + gp_surf(3, i) = 0._wp + if (reacting) gp_surf(3, i) = mdot_s + + end subroutine s_record_gp_surface + + !> Axis (1, 2 or 3) of a cylinder IB: the one length that is set. + pure integer function f_cylinder_axis(ib_patch) + + $:GPU_ROUTINE(parallelism='[seq]') + + type(ib_patch_parameters), intent(in) :: ib_patch + + f_cylinder_axis = 3 + if (ib_patch%length_x > 0._wp) f_cylinder_axis = 1 + if (ib_patch%length_y > 0._wp) f_cylinder_axis = 2 + + end function f_cylinder_axis + + !> Write this rank's surface records to D/ib_surface__.dat, one line per weighted ghost point. + impure subroutine s_write_ib_surface(save_count) + + integer, intent(in) :: save_count + character(LEN=path_len + 2*name_len) :: file_loc + real(wp) :: x(3) + integer :: i, unit + + if (.not. allocated(gp_surf)) return + + if (num_gps > 0) then + $:GPU_UPDATE(host='[gp_surf(:, 1:num_gps), ghost_points(1:num_gps)]') + end if + + write (file_loc, '(A,I0,A,I0,A)') trim(case_dir) // '/D/ib_surface_', proc_rank, '_', save_count, '.dat' + open (newunit=unit, file=trim(file_loc), status='replace', action='write') + do i = 1, num_gps + if (.not. gp_surf(1, i) > 0._wp) cycle + x = 0._wp + x(1) = x_cc(ghost_points(i)%loc(1)) + x(2) = y_cc(ghost_points(i)%loc(2)) + if (num_dims == 3) x(3) = z_cc(ghost_points(i)%loc(3)) + write (unit, '(3ES16.8,1X,I0,3ES15.6,3ES16.8)') x, patch_ib(ghost_points(i)%ib_patch_id)%gbl_patch_id, & + & ghost_points(i)%levelset_norm, gp_surf(:,i) + end do + close (unit) + + end subroutine s_write_ib_surface + !> Compute the image points for each ghost point impure subroutine s_compute_image_points() @@ -2020,6 +2114,9 @@ contains if (allocated(ghost_points)) then @:DEALLOCATE(ghost_points) end if + if (allocated(gp_surf)) then + @:DEALLOCATE(gp_surf) + end if if (collision_model > 0) call s_finalize_collisions_module() #ifdef MFC_MPI if (allocated(ib_nbrs)) deallocate (ib_nbrs, ib_send_counts, ib_reqs, ib_stats, ib_send_bufs, ib_recv_bufs) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 4933921df..7b0d320c7 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -807,6 +807,7 @@ contains ! Write IB kinematic state for restart if (ib) call s_write_ib_state_file(save_count) + if (ib .and. ib_surface_wrt) call s_write_ib_surface(save_count) call nvtxEndRange call cpu_time(finish) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 36fce10a0..41f497676 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -799,6 +799,8 @@ def check_ibm(self): self.prohibit(ib_state_wrt and not ib, "ib_state_wrt requires ib to be enabled") ib_force_wrt = self.get("ib_force_wrt", False) self.prohibit(ib_force_wrt and not ib, "ib_force_wrt requires ib to be enabled") + ib_surface_wrt = self.get("ib_surface_wrt", "F") == "T" + self.prohibit(ib_surface_wrt and not (ib and self.get("chemistry", "F") == "T"), "ib_surface_wrt requires ib and chemistry") ib_force_stride = self.get("ib_force_stride", 1) self.prohibit(ib_force_stride < 1, "ib_force_stride must be >= 1") diff --git a/toolchain/mfc/lint_docs.py b/toolchain/mfc/lint_docs.py index e86b24d1d..f772ffa6a 100644 --- a/toolchain/mfc/lint_docs.py +++ b/toolchain/mfc/lint_docs.py @@ -72,6 +72,9 @@ "m_constants", # Build/run target name (not a case param) "pre_process", + # Output file column names (ib_surface_wrt) + "area", + "mdot", } # Docs to check for parameter references, with per-file skip sets diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 159dd9aae..d144181c3 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -703,7 +703,7 @@ def _load(): _r("precision", INT, {"output"}) _r("format", INT, {"output"}) _r("ib_force_stride", INT, {"output", "ib"}) - for n in ["parallel_io", "file_per_process", "run_time_info", "prim_vars_wrt", "cons_vars_wrt", "fft_wrt", "ib_state_wrt", "ib_force_wrt"]: + for n in ["parallel_io", "file_per_process", "run_time_info", "prim_vars_wrt", "cons_vars_wrt", "fft_wrt", "ib_state_wrt", "ib_force_wrt", "ib_surface_wrt"]: _r(n, LOG, {"output"}) for n in [ "schlieren_wrt", @@ -1356,6 +1356,7 @@ def _decl(targets: set, *names: str) -> None: "ib_state_wrt", "ib_force_wrt", "ib_force_stride", + "ib_surface_wrt", "avg_state", "alt_soundspeed", "mixture_err", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index 7dfb46f3e..4ef88b641 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -137,6 +137,7 @@ # Immersed boundaries "ib": "Enable immersed boundary method", "ib_force_wrt": "Record the immersed-boundary force history to D/ib_forces.dat (default off)", + "ib_surface_wrt": "Write per-surface-point wall temperature and gasified mass flux of thermal/reacting IBs to D/ib_surface__.dat (default off)", "ib_force_stride": "Write the per-step immersed-boundary force record every N steps (default 1)", "num_ibs": "Number of immersed boundary patches", "num_stl_models": "Number of STL/OBJ model entries in the stl_models array", diff --git a/toolchain/mfc/test_case_validator.py b/toolchain/mfc/test_case_validator.py index c0170452d..4ba04c6a8 100644 --- a/toolchain/mfc/test_case_validator.py +++ b/toolchain/mfc/test_case_validator.py @@ -625,6 +625,10 @@ def test_twall_must_lie_in_the_tabulated_range(self): self.assertRejects({**self.case(thermal_bc=1, Twall=9000.0), "chemistry": "T"}, "Twall must be within") self.assertAccepts({**self.case(thermal_bc=1, Twall=210.0), "chemistry": "T"}) + def test_surface_output_needs_chemistry(self): + self.assertRejects({**self.case(thermal_bc=1, Twall=1200.0), "ib_surface_wrt": "T"}, "ib_surface_wrt requires ib and chemistry") + self.assertAccepts({**self.case(thermal_bc=1, Twall=1200.0), "chemistry": "T", "ib_surface_wrt": "T"}) + def test_the_energy_balance_needs_a_reacting_surface(self): self.assertRejects({**self.case(thermal_bc=2), "chemistry": "T"}, "thermal_bc = 2 requires surface_reaction = 1") self.assertAccepts({**self.case(thermal_bc=2, surface_reaction=1), "chemistry": "T", "surface_cantera_file": "s.yaml", "surface_phase": "surf"})