!> @file mesh_1d.f90 !> @brief 1D grid metric layer: node coordinates and the dx/dξ Jacobian. !! !! The 1D solver is node-based finite-difference WENO. On a non-uniform grid the !! physical coordinate is mapped to a uniform computational grid ξ (Δξ = 1); the !! residual divides the flux difference by the per-node Jacobian J_i = dx/dξ !! rather than a scalar dx. This module owns the node coordinates and J. !! !! Dependency rule: this module must NOT `use config` (config uses this module). module mesh_1d use precision, only: wp use domain_decomposition, only: decomp_t implicit none private public :: mesh_1d_t public :: read_node_coords public :: build_mesh_global public :: build_mesh_local public :: build_mesh_cellcentered_1d public :: build_mesh_cellcentered_global public :: build_mesh_cellcentered_local !> Grid metrics for the 1D node-based scheme. !! !! The nodal fields (`x_global`, `jac_global`, `x_node`, `jac`) are used by !! the FDM path. The cell-centered fields (`x_cell`, `dx_cell`) are used by !! the FVM path and are allocated by `build_mesh_cellcentered_1d`; they are !! unallocated (and unused) on the FDM path. type :: mesh_1d_t logical :: uniform = .true. integer :: n_global = 0 real(wp) :: dx_uniform = 0.0_wp !< constant spacing (uniform path only) real(wp) :: h_min = 0.0_wp !< min node spacing (CFL) real(wp) :: h_left = 0.0_wp !< first interval width (Neumann-gradient BC) real(wp) :: h_right = 0.0_wp !< last interval width (Neumann-gradient BC) real(wp), allocatable :: x_global(:) !< global node coords (1:n_global) real(wp), allocatable :: jac_global(:) !< global dx/dξ (1:n_global) real(wp), allocatable :: x_node(:) !< local node coords (1:n_local) real(wp), allocatable :: jac(:) !< local dx/dξ (1:n_local) !> Cell-center coordinates for the FVM path. !! Allocated with bounds (1-h:n_cell+h) by build_mesh_cellcentered_1d. real(wp), allocatable :: x_cell(:) !> Cell widths for the FVM path (uniform: all equal dx). !! Allocated with bounds (1-h:n_cell+h) by build_mesh_cellcentered_1d. real(wp), allocatable :: dx_cell(:) !> Replicated GLOBAL cell-center coordinates for a NON-UNIFORM FVM grid !! (1:n_cell). Populated by build_mesh_cellcentered_global from the face !! positions stored in x_global; sliced per rank by !! build_mesh_cellcentered_local. Unallocated on the uniform FVM path, !! which recomputes cell centers analytically from x_left/x_right. real(wp), allocatable :: x_cell_global(:) !> Replicated GLOBAL cell widths for a NON-UNIFORM FVM grid (1:n_cell), !! the cell-centered twin of jac_global. See x_cell_global. real(wp), allocatable :: dx_cell_global(:) end type mesh_1d_t contains !> Read strictly-increasing node coordinates from a text file. !! One value per line; blank lines and lines beginning with '#' are ignored. subroutine read_node_coords(path, coords, n, is_ok, message) character(len=*), intent(in) :: path real(wp), allocatable, intent(out) :: coords(:) integer, intent(out) :: n logical, intent(out) :: is_ok character(len=*), intent(out) :: message integer :: u, info, nline, k, dstat real(wp) :: val, prev character(len=512) :: line is_ok = .false. message = '' n = 0 prev = 0.0_wp open (newunit=u, file=trim(path), status='old', action='read', iostat=info) if (info /= 0) then message = 'mesh_1d: cannot open grid_file "'//trim(path)//'"' return end if ! Pass 1: count data lines. nline = 0 do read (u, '(A)', iostat=info) line if (info /= 0) exit line = adjustl(line) if (len_trim(line) == 0) cycle if (line(1:1) == '#') cycle nline = nline + 1 end do if (nline < 2) then close (u, iostat=info) message = 'mesh_1d: grid_file must contain at least 2 node coordinates' return end if allocate (coords(nline), stat=info) if (info /= 0) error stop 'mesh_1d: coords allocation failed' ! Pass 2: parse and check strict monotonicity. rewind (u) k = 0 do read (u, '(A)', iostat=info) line if (info /= 0) exit line = adjustl(line) if (len_trim(line) == 0) cycle if (line(1:1) == '#') cycle read (line, *, iostat=info) val if (info /= 0) then close (u, iostat=info) deallocate (coords, stat=dstat) message = 'mesh_1d: malformed number in grid_file' if (dstat /= 0) message = trim(message)//' (note: coords deallocate failed)' return end if k = k + 1 if (k > 1 .and. val <= prev) then close (u, iostat=info) deallocate (coords, stat=dstat) write (message, '(A,I0)') 'mesh_1d: grid_file not strictly increasing at node ', k if (dstat /= 0) message = trim(message)//' (note: coords deallocate failed)' return end if coords(k) = val prev = val end do close (u, iostat=info) if (k /= nline) then deallocate (coords, stat=dstat) message = 'mesh_1d: grid_file changed during read (line-count mismatch)' if (dstat /= 0) message = trim(message)//' (note: coords deallocate failed)' return end if n = k is_ok = .true. end subroutine read_node_coords !> Build the replicated global mesh; for file grids, derive n_cell/endpoints. subroutine build_mesh_global(mesh, grid_type, grid_file, n_cell, x_left, x_right, is_ok, message) type(mesh_1d_t), intent(inout) :: mesh character(len=*), intent(in) :: grid_type, grid_file integer, intent(inout) :: n_cell real(wp), intent(inout) :: x_left, x_right logical, intent(out) :: is_ok character(len=*), intent(out) :: message integer :: k, n, astat real(wp), allocatable :: coords(:) is_ok = .true. message = '' if (trim(grid_type) == 'file') then call read_node_coords(grid_file, coords, n, is_ok, message) if (.not. is_ok) return mesh % uniform = .false. mesh % n_global = n mesh % dx_uniform = 0.0_wp ! no single spacing on a non-uniform grid call move_alloc(coords, mesh % x_global) n_cell = n - 1 ! file is authoritative x_left = mesh % x_global(1) x_right = mesh % x_global(n) else if (trim(grid_type) == 'uniform') then if (n_cell <= 0) then is_ok = .false. message = 'mesh_1d: uniform grid requires n_cell > 0' return end if mesh % uniform = .true. n = n_cell + 1 mesh % n_global = n mesh % dx_uniform = (x_right - x_left) / real(n_cell, wp) if (allocated(mesh % x_global)) then deallocate (mesh % x_global, stat=astat) if (astat /= 0) error stop 'mesh_1d: x_global deallocate failed' end if allocate (mesh % x_global(n), stat=astat) if (astat /= 0) error stop 'mesh_1d: x_global allocation failed' do k = 1, n mesh % x_global(k) = x_left + mesh % dx_uniform * real(k - 1, wp) end do else ! Reject unknown grid_type rather than silently treating it as uniform. ! validate_config also guards this, but build_mesh_global is called ! directly by tests that bypass validation. is_ok = .false. message = "mesh_1d: unknown grid_type '"//trim(grid_type)// & "' (expected 'uniform' or 'file')" return end if ! Jacobian J_k = dx/dξ on the uniform ξ grid (Δξ = 1). if (allocated(mesh % jac_global)) then deallocate (mesh % jac_global, stat=astat) if (astat /= 0) error stop 'mesh_1d: jac_global deallocate failed' end if allocate (mesh % jac_global(n), stat=astat) if (astat /= 0) error stop 'mesh_1d: jac_global allocation failed' if (mesh % uniform) then mesh % jac_global = mesh % dx_uniform ! exact -> bit-for-bit with /dx mesh % h_min = mesh % dx_uniform mesh % h_left = mesh % dx_uniform mesh % h_right = mesh % dx_uniform else mesh % jac_global(1) = mesh % x_global(2) - mesh % x_global(1) mesh % jac_global(n) = mesh % x_global(n) - mesh % x_global(n - 1) do k = 2, n - 1 mesh % jac_global(k) = 0.5_wp * (mesh % x_global(k + 1) - mesh % x_global(k - 1)) end do mesh % h_left = mesh % x_global(2) - mesh % x_global(1) mesh % h_right = mesh % x_global(n) - mesh % x_global(n - 1) mesh % h_min = mesh % h_left do k = 2, n - 1 mesh % h_min = min(mesh % h_min, mesh % x_global(k + 1) - mesh % x_global(k)) end do end if end subroutine build_mesh_global !> Slice the replicated global mesh into this rank's local node arrays. subroutine build_mesh_local(mesh, decomp) type(mesh_1d_t), intent(inout) :: mesh type(decomp_t), intent(in) :: decomp integer :: ipt, g, astat if (.not. allocated(mesh % x_global) .or. .not. allocated(mesh % jac_global)) & error stop 'mesh_1d: build_mesh_local called before build_mesh_global' if (allocated(mesh % x_node)) then deallocate (mesh % x_node, stat=astat) if (astat /= 0) error stop 'mesh_1d: x_node deallocate failed' end if if (allocated(mesh % jac)) then deallocate (mesh % jac, stat=astat) if (astat /= 0) error stop 'mesh_1d: jac deallocate failed' end if allocate (mesh % x_node(decomp % n_local), stat=astat) if (astat /= 0) error stop 'mesh_1d: x_node allocation failed' allocate (mesh % jac(decomp % n_local), stat=astat) if (astat /= 0) error stop 'mesh_1d: jac allocation failed' do ipt = 1, decomp % n_local g = decomp % i_first_global - 1 + ipt mesh % x_node(ipt) = mesh % x_global(g) mesh % jac(ipt) = mesh % jac_global(g) end do end subroutine build_mesh_local !> Build cell-centered coordinates and widths for the FVM path. !! !! For a uniform domain `[x_left, x_right]` divided into `n_cell` cells: !! - `dx = (x_right - x_left) / n_cell` !! - interior cell centers: `x_cell(i) = x_left + (i - 0.5) * dx` (i = 1..n_cell) !! - ghost centers extrapolated by `dx`: `x_cell(i) = x_left + (i - 0.5) * dx` !! (the same formula applies for halo indices 1-h..0 and n_cell+1..n_cell+h) !! - `dx_cell(i) = dx` for all i in `1-halo_width:n_cell+halo_width` !! !! Both arrays are allocated with lower bound `1 - halo_width`. !! !! @param mesh Mesh object to populate (x_cell and dx_cell are set). !! @param n_cell Number of interior cells. !! @param x_left Left domain boundary coordinate. !! @param x_right Right domain boundary coordinate. !! @param halo_width Number of ghost-cell layers on each side (h >= 1). subroutine build_mesh_cellcentered_1d(mesh, n_cell, x_left, x_right, halo_width) type(mesh_1d_t), intent(inout) :: mesh integer, intent(in) :: n_cell real(wp), intent(in) :: x_left, x_right integer, intent(in) :: halo_width integer :: i, lo, hi, astat real(wp) :: dx lo = 1 - halo_width hi = n_cell + halo_width dx = (x_right - x_left) / real(n_cell, wp) if (allocated(mesh % x_cell)) then deallocate (mesh % x_cell, stat=astat) if (astat /= 0) error stop 'mesh_1d: x_cell deallocate failed' end if if (allocated(mesh % dx_cell)) then deallocate (mesh % dx_cell, stat=astat) if (astat /= 0) error stop 'mesh_1d: dx_cell deallocate failed' end if allocate (mesh % x_cell(lo:hi), stat=astat) if (astat /= 0) error stop 'mesh_1d: x_cell allocation failed' allocate (mesh % dx_cell(lo:hi), stat=astat) if (astat /= 0) error stop 'mesh_1d: dx_cell allocation failed' do i = lo, hi mesh % x_cell(i) = x_left + (real(i, wp) - 0.5_wp) * dx mesh % dx_cell(i) = dx end do end subroutine build_mesh_cellcentered_1d !> Build the replicated GLOBAL cell-centered mesh for a NON-UNIFORM FVM grid. !! !! The cell-centered twin of the nodal arrays built by build_mesh_global: where !! that routine fills x_global/jac_global, this fills x_cell_global and !! dx_cell_global. For an FVM file grid the grid_file nodes ARE the cell !! faces, and build_mesh_global has already loaded those n_cell+1 faces into !! `mesh % x_global`. This routine derives, for each interior cell !! i = 1..n_cell: !! - center x_cell_global(i) = 0.5 * (face(i) + face(i+1)) !! - width dx_cell_global(i) = face(i+1) - face(i) !! and is the authoritative source of non-uniform cell geometry, later sliced !! per rank by build_mesh_cellcentered_local. !! !! It is intentionally NOT called on the uniform FVM path: there the local !! builder recomputes cell centers/widths analytically from x_left/x_right so !! the uniform result stays bit-for-bit identical to the original formula. !! !! @param mesh Mesh object; x_global (faces) must already be populated. !! @param n_cell Number of interior cells (= number of faces - 1). subroutine build_mesh_cellcentered_global(mesh, n_cell) type(mesh_1d_t), intent(inout) :: mesh integer, intent(in) :: n_cell integer :: i, astat if (.not. allocated(mesh % x_global)) & error stop 'mesh_1d: build_mesh_cellcentered_global requires x_global (cell faces)' if (size(mesh % x_global) < n_cell + 1) & error stop 'mesh_1d: build_mesh_cellcentered_global: x_global has fewer than n_cell+1 faces' if (allocated(mesh % x_cell_global)) then deallocate (mesh % x_cell_global, stat=astat) if (astat /= 0) error stop 'mesh_1d: x_cell_global deallocate failed' end if if (allocated(mesh % dx_cell_global)) then deallocate (mesh % dx_cell_global, stat=astat) if (astat /= 0) error stop 'mesh_1d: dx_cell_global deallocate failed' end if allocate (mesh % x_cell_global(n_cell), stat=astat) if (astat /= 0) error stop 'mesh_1d: x_cell_global allocation failed' allocate (mesh % dx_cell_global(n_cell), stat=astat) if (astat /= 0) error stop 'mesh_1d: dx_cell_global allocation failed' do i = 1, n_cell mesh % x_cell_global(i) = 0.5_wp * (mesh % x_global(i) + mesh % x_global(i + 1)) mesh % dx_cell_global(i) = mesh % x_global(i + 1) - mesh % x_global(i) end do end subroutine build_mesh_cellcentered_global !> Build this rank's local cell-centered coordinates and widths (FVM path). !! !! The MPI twin of `build_mesh_cellcentered_1d`: where that routine builds the !! replicated GLOBAL cell-centered mesh, this slices out the per-rank LOCAL !! cell-centered mesh from the decomposition, mirroring how `build_mesh_local` !! slices the nodal mesh. For local cell `i` (1..n_local) the 1-based global !! cell index is `decomp % i_first_global - 1 + i`; the uniform cell width is !! `dx = (x_right - x_left) / decomp % n_global` (the FVM decomposition counts !! CELLS, so `n_global` is the global cell count); the cell center is !! `x_left + (global_i - 0.5) * dx`. !! !! UNIFORM grids: the ghost layers `1-h .. 0` and `n_local+1 .. n_local+h` are !! filled with the same center formula evaluated at their (possibly !! out-of-domain) global cell index, so the global-edge ghosts extrapolate !! linearly by `dx`. At np=1 (`i_first_global == 1`, `n_local == n_global`) !! this reduces exactly to !! `build_mesh_cellcentered_1d(mesh, n_global, x_left, x_right, halo_width)`. !! !! NON-UNIFORM grids (`.not. mesh % uniform`, i.e. an FVM file grid): interior !! cells are sliced from the replicated global arrays built by !! build_mesh_cellcentered_global — `x_cell(i) = x_cell_global(g)`, !! `dx_cell(i) = dx_cell_global(g)` with `g = i_first_global - 1 + i`, mirroring !! build_mesh_local's nodal slice. Ghost cells beyond the global domain !! (`g < 1` or `g > n_global`) are filled by linear extrapolation of the edge !! face positions: the nearest edge cell width is held constant and the center !! steps out by that width. These ghost geometry values are NOT used by the !! residual (which only divides interior cells `1..n_local` by their own !! dx_cell) nor by the uniform-coefficient reconstruction (which is !! geometry-agnostic); they are filled only so the arrays are consistent. !! !! @param mesh Mesh object to populate (x_cell and dx_cell are set). !! @param decomp This rank's decomposition (cell-counted for the FVM path). !! @param x_left Left domain boundary coordinate. !! @param x_right Right domain boundary coordinate. subroutine build_mesh_cellcentered_local(mesh, decomp, x_left, x_right) type(mesh_1d_t), intent(inout) :: mesh type(decomp_t), intent(in) :: decomp real(wp), intent(in) :: x_left, x_right integer :: i, g, lo, hi, h, ng, astat real(wp) :: dx h = decomp % halo_width lo = 1 - h hi = decomp % n_local + h ng = decomp % n_global if (allocated(mesh % x_cell)) then deallocate (mesh % x_cell, stat=astat) if (astat /= 0) error stop 'mesh_1d: x_cell deallocate failed' end if if (allocated(mesh % dx_cell)) then deallocate (mesh % dx_cell, stat=astat) if (astat /= 0) error stop 'mesh_1d: dx_cell deallocate failed' end if allocate (mesh % x_cell(lo:hi), stat=astat) if (astat /= 0) error stop 'mesh_1d: x_cell allocation failed' allocate (mesh % dx_cell(lo:hi), stat=astat) if (astat /= 0) error stop 'mesh_1d: dx_cell allocation failed' if (mesh % uniform) then dx = (x_right - x_left) / real(ng, wp) do i = lo, hi g = decomp % i_first_global - 1 + i mesh % x_cell(i) = x_left + (real(g, wp) - 0.5_wp) * dx mesh % dx_cell(i) = dx end do else if (.not. allocated(mesh % x_cell_global) .or. .not. allocated(mesh % dx_cell_global)) & error stop 'mesh_1d: non-uniform build_mesh_cellcentered_local requires the global cell mesh' do i = lo, hi g = decomp % i_first_global - 1 + i if (g >= 1 .and. g <= ng) then mesh % x_cell(i) = mesh % x_cell_global(g) mesh % dx_cell(i) = mesh % dx_cell_global(g) else if (g < 1) then ! Left ghost: extend beyond the first cell by its width. mesh % dx_cell(i) = mesh % dx_cell_global(1) mesh % x_cell(i) = mesh % x_cell_global(1) - real(1 - g, wp) * mesh % dx_cell_global(1) else ! Right ghost: extend beyond the last cell by its width. mesh % dx_cell(i) = mesh % dx_cell_global(ng) mesh % x_cell(i) = mesh % x_cell_global(ng) + real(g - ng, wp) * mesh % dx_cell_global(ng) end if end do end if end subroutine build_mesh_cellcentered_local end module mesh_1d