!> @file mesh_2d.f90 !> @brief 2D single-block structured (curvilinear) grid metric layer. !! !! Owns node coordinates and the derived finite-volume metrics for the 2D Euler !! solver: cell areas (vol), face area-normal vectors (s_xi, s_eta), and cell !! centres (xc, yc). Geometry is read from a single-block Plot3D grid file !! (ASCII or binary). On a uniform Cartesian grid the metrics are unused — the !! residual keeps its scalar dx/dy fast path (bit-for-bit). Phase A builds the !! whole-domain (serial) metrics; local slicing for MPI is added in Phase D. !! !! Dependency rule: this module must NOT `use config_2d` (config_2d uses it). module mesh_2d use precision, only: wp use, intrinsic :: iso_fortran_env, only: int64 use, intrinsic :: ieee_arithmetic, only: ieee_is_finite implicit none private public :: mesh_2d_t public :: read_plot3d_2d public :: check_plot3d_dims public :: build_mesh_2d_global public :: build_mesh_2d_uniform_nodal public :: compute_metrics_2d public :: compute_metrics_fdm_2d public :: build_mesh_2d_fdm_curvilinear public :: build_mesh_2d_local public :: MESH2D_MAX_DIM !> Grid metrics for the 2D structured scheme. type :: mesh_2d_t logical :: uniform = .true. !< true => scalar dx/dy fast path integer :: nx = 0, ny = 0 !< cell count per axis (FVM) or node count (FDM nodal mesh) real(wp) :: dx_uniform = 0.0_wp, dy_uniform = 0.0_wp real(wp), allocatable :: x_node(:, :) !< node coords; shape (nx+1,ny+1) for FVM; (mesh%nx,mesh%ny) for the FDM nodal mesh real(wp), allocatable :: y_node(:, :) real(wp), allocatable :: vol(:, :) !< cell area V_{i,j} (nx, ny) real(wp), allocatable :: s_xi(:, :, :) !< xi-face vectors (2, nx+1, ny) real(wp), allocatable :: s_eta(:, :, :)!< eta-face vectors (2, nx, ny+1) real(wp), allocatable :: xc(:, :) !< cell-centre x (nx, ny) real(wp), allocatable :: yc(:, :) !< cell-centre y (nx, ny) !> Nodal FDM transformation metrics (allocated only for curvilinear FDM). real(wp), allocatable :: sx_xi(:, :, :) !< S_xi=(y_eta,-x_eta) at xi-faces (2, nx+1, ny) real(wp), allocatable :: sx_eta(:, :, :) !< S_eta=(-y_xi,x_xi) at eta-faces (2, nx, ny+1) real(wp), allocatable :: jac(:, :) !< Jacobian J>0 at nodes (nx, ny) logical :: fdm_curvilinear = .false. !< true when nodal metrics are present end type mesh_2d_t !> Sanity caps on Plot3D header dimensions. Per-axis node cap mirrors the !! spirit of the 1D n_cell schema bound; the total-node cap is checked via an !! int64 product so a pathological header cannot overflow the default-integer !! ni*nj before the allocation's stat= is consulted. integer, parameter :: MESH2D_MAX_DIM = 1000000 integer(int64), parameter :: MESH2D_MAX_POINTS = 100000000_int64 contains !> Read a single-block 2D Plot3D grid. Header `ni nj`, then all x(i,j) !! (i fastest), then all y(i,j). ASCII (binary=.false.) or Fortran !! unformatted sequential (binary=.true.). subroutine read_plot3d_2d(path, binary, xg, yg, ni, nj, is_ok, message) character(len=*), intent(in) :: path logical, intent(in) :: binary real(wp), allocatable, intent(out) :: xg(:, :), yg(:, :) integer, intent(out) :: ni, nj logical, intent(out) :: is_ok character(len=*), intent(out) :: message integer :: u, info, astat, i, j is_ok = .false. message = '' ni = 0; nj = 0 if (binary) then open (newunit=u, file=trim(path), status='old', action='read', & form='unformatted', access='sequential', iostat=info) if (info /= 0) then message = 'mesh_2d: cannot open binary grid_file "'//trim(path)//'"' return end if read (u, iostat=info) ni, nj if (info /= 0) then close (u, iostat=info) message = 'mesh_2d: malformed Plot3D binary header' return end if block logical :: dims_ok call check_plot3d_dims(ni, nj, dims_ok, message) if (.not. dims_ok) then close (u, iostat=info) return end if end block allocate (xg(ni, nj), yg(ni, nj), stat=astat) if (astat /= 0) then close (u, iostat=info) message = 'mesh_2d: allocation failed for grid arrays' return end if read (u, iostat=info) ((xg(i, j), i=1, ni), j=1, nj) if (info == 0) read (u, iostat=info) ((yg(i, j), i=1, ni), j=1, nj) close (u, iostat=astat) if (info /= 0) then deallocate (xg, yg, stat=astat) message = 'mesh_2d: malformed/short Plot3D binary coordinate block'// & ' (also check host endianness for binary grids)' return end if if (.not. (all_finite_2d(xg) .and. all_finite_2d(yg))) then deallocate (xg, yg, stat=astat) message = 'mesh_2d: non-finite (NaN/Inf) coordinate in binary Plot3D grid'// & ' (truncated record or wrong endianness?)' return end if is_ok = .true. return end if ! ASCII path. open (newunit=u, file=trim(path), status='old', action='read', iostat=info) if (info /= 0) then message = 'mesh_2d: cannot open grid_file "'//trim(path)//'"' return end if read (u, *, iostat=info) ni, nj if (info /= 0) then close (u, iostat=info) message = 'mesh_2d: malformed Plot3D header (expected "ni nj")' return end if block logical :: dims_ok call check_plot3d_dims(ni, nj, dims_ok, message) if (.not. dims_ok) then close (u, iostat=info) return end if end block allocate (xg(ni, nj), yg(ni, nj), stat=astat) if (astat /= 0) then close (u, iostat=info) message = 'mesh_2d: allocation failed for grid arrays' return end if read (u, *, iostat=info) ((xg(i, j), i=1, ni), j=1, nj) if (info == 0) read (u, *, iostat=info) ((yg(i, j), i=1, ni), j=1, nj) close (u, iostat=astat) if (info /= 0) then deallocate (xg, yg, stat=astat) message = 'mesh_2d: malformed/short Plot3D coordinate block' return end if if (.not. (all_finite_2d(xg) .and. all_finite_2d(yg))) then deallocate (xg, yg, stat=astat) message = 'mesh_2d: non-finite (NaN/Inf) coordinate in Plot3D grid' return end if is_ok = .true. end subroutine read_plot3d_2d !> True iff every element of a 2D coordinate array is finite (no NaN/Inf). pure logical function all_finite_2d(a) result(ok) real(wp), intent(in) :: a(:, :) ok = all(ieee_is_finite(a)) end function all_finite_2d !> Validate a Plot3D header: lower bound (>= 2 nodes per axis) and upper !! sanity caps (per-axis and total node count). Returns is_ok=.false. with a !! descriptive message on violation. subroutine check_plot3d_dims(ni, nj, is_ok, message) integer, intent(in) :: ni, nj logical, intent(out) :: is_ok character(len=*), intent(out) :: message is_ok = .false. message = '' ! Per-axis and total-node caps are checked in separate branches so the ! diagnostic names the cap that was actually exceeded. (The int64 product ! cannot overflow: both operands are promoted to int64 before multiplying, ! and any int32 ni*nj is <= ~4.6e18 < huge(int64).) if (ni < 2 .or. nj < 2) then write (message, '(A,I0,A,I0)') 'mesh_2d: grid too small ni=', ni, ' nj=', nj else if (ni > MESH2D_MAX_DIM .or. nj > MESH2D_MAX_DIM) then write (message, '(A,I0,A,I0,A,I0,A)') 'mesh_2d: grid dimension exceeds per-axis cap ni=', & ni, ' nj=', nj, ' (max ', MESH2D_MAX_DIM, ' per axis)' else if (int(ni, int64) * int(nj, int64) > MESH2D_MAX_POINTS) then write (message, '(A,I0,A,I0,A,I0,A)') 'mesh_2d: grid total node count exceeds cap ni=', & ni, ' nj=', nj, ' (max ', MESH2D_MAX_POINTS, ' total)' else is_ok = .true. end if end subroutine check_plot3d_dims !> Build the global (whole-domain) mesh from either a uniform specification or a !! Plot3D grid file. For 'uniform', nx/ny are read from the caller. For !! 'plot3d', ni/nj are read from the file and nx/ny are written back !! (intent(inout)) so the caller knows the derived sizes. Calls !! compute_metrics_2d at the end and propagates its status — a metric failure !! (e.g. a folded or zero-area cell) fails the build. subroutine build_mesh_2d_global(mesh, grid_type, grid_file, grid_binary, & nx, ny, x_left, x_right, y_left, y_right, & is_ok, message) type(mesh_2d_t), intent(inout) :: mesh character(len=*), intent(in) :: grid_type, grid_file logical, intent(in) :: grid_binary integer, intent(inout) :: nx, ny real(wp), intent(in) :: x_left, x_right, y_left, y_right logical, intent(out) :: is_ok character(len=*), intent(out) :: message integer :: ni, nj, i, j, astat real(wp), allocatable :: xg(:, :), yg(:, :) is_ok = .true. message = '' if (trim(grid_type) == 'plot3d') then call read_plot3d_2d(grid_file, grid_binary, xg, yg, ni, nj, is_ok, message) if (.not. is_ok) return mesh % uniform = .false. mesh % nx = ni - 1 mesh % ny = nj - 1 mesh % dx_uniform = 0.0_wp mesh % dy_uniform = 0.0_wp call move_alloc(xg, mesh % x_node) call move_alloc(yg, mesh % y_node) nx = mesh % nx ! file is authoritative ny = mesh % ny else if (trim(grid_type) == 'uniform') then if (nx <= 0 .or. ny <= 0) then is_ok = .false. message = 'mesh_2d: uniform grid requires nx > 0 and ny > 0' return end if mesh % uniform = .true. mesh % nx = nx mesh % ny = ny mesh % dx_uniform = (x_right - x_left) / real(nx, wp) mesh % dy_uniform = (y_right - y_left) / real(ny, wp) ni = nx + 1; nj = ny + 1 if (allocated(mesh % x_node)) then deallocate (mesh % x_node, stat=astat) if (astat /= 0) error stop 'mesh_2d: x_node deallocate failed' end if if (allocated(mesh % y_node)) then deallocate (mesh % y_node, stat=astat) if (astat /= 0) error stop 'mesh_2d: y_node deallocate failed' end if allocate (mesh % x_node(ni, nj), mesh % y_node(ni, nj), stat=astat) if (astat /= 0) then is_ok = .false. message = 'mesh_2d: allocation failed for uniform node arrays' return end if do j = 1, nj do i = 1, ni mesh % x_node(i, j) = x_left + mesh % dx_uniform * real(i - 1, wp) mesh % y_node(i, j) = y_left + mesh % dy_uniform * real(j - 1, wp) end do end do else is_ok = .false. message = "mesh_2d: unknown grid_type '"//trim(grid_type)// & "' (expected 'uniform' or 'plot3d')" return end if call compute_metrics_2d(mesh, is_ok, message) end subroutine build_mesh_2d_global !> Compute finite-volume metrics from node coordinates already in mesh%x_node / !! mesh%y_node. Allocates (with stat=) and fills vol(nx,ny), s_xi(2,nx+1,ny), !! s_eta(2,nx,ny+1), xc(nx,ny), yc(nx,ny). Returns is_ok=.false. and a !! diagnostic message naming the first (i,j) with non-positive cell area. !! !! Metric definitions !! ------------------ !! * xi-face vector at column iface, row j: edge (iface,j)->(iface,j+1). !! s_xi = (dy, -dx), i.e. (+dy, -dx) where dx/dy are the edge increments. !! * eta-face vector at row jface, column i: edge (i,jface)->(i+1,jface). !! s_eta = (-dy, dx). !! * Cell area via diagonal cross-product (exact for arbitrary quads): !! vol = 0.5 * [(d1x*d2y) - (d2x*d1y)] where d1,d2 are the two diagonals. !! * Cell centre = arithmetic mean of four corners. subroutine compute_metrics_2d(mesh, is_ok, message) type(mesh_2d_t), intent(inout) :: mesh logical, intent(out) :: is_ok character(len=*), intent(out) :: message integer :: nx, ny, i, j, astat real(wp) :: dx, dy is_ok = .true. message = '' nx = mesh % nx; ny = mesh % ny if (allocated(mesh % vol)) then deallocate (mesh % vol, stat=astat) if (astat /= 0) error stop 'mesh_2d: vol deallocate failed' end if if (allocated(mesh % s_xi)) then deallocate (mesh % s_xi, stat=astat) if (astat /= 0) error stop 'mesh_2d: s_xi deallocate failed' end if if (allocated(mesh % s_eta)) then deallocate (mesh % s_eta, stat=astat) if (astat /= 0) error stop 'mesh_2d: s_eta deallocate failed' end if if (allocated(mesh % xc)) then deallocate (mesh % xc, stat=astat) if (astat /= 0) error stop 'mesh_2d: xc deallocate failed' end if if (allocated(mesh % yc)) then deallocate (mesh % yc, stat=astat) if (astat /= 0) error stop 'mesh_2d: yc deallocate failed' end if allocate (mesh % vol(nx, ny), mesh % xc(nx, ny), mesh % yc(nx, ny), & mesh % s_xi(2, nx + 1, ny), mesh % s_eta(2, nx, ny + 1), stat=astat) if (astat /= 0) then is_ok = .false. message = 'mesh_2d: allocation failed for metric arrays' return end if ! xi-face vectors at node columns iface=1..nx+1, edge (iface,j)->(iface,j+1). do j = 1, ny do i = 1, nx + 1 dx = mesh % x_node(i, j + 1) - mesh % x_node(i, j) dy = mesh % y_node(i, j + 1) - mesh % y_node(i, j) mesh % s_xi(1, i, j) = dy mesh % s_xi(2, i, j) = -dx if (mesh % s_xi(1, i, j)**2 + mesh % s_xi(2, i, j)**2 <= 0.0_wp) then is_ok = .false. write (message, '(A,I0,A,I0,A)') 'mesh_2d: zero-length xi-face at (', i, ',', j, ')' return end if end do end do ! eta-face vectors at node rows jface=1..ny+1, edge (i,jface)->(i+1,jface). do j = 1, ny + 1 do i = 1, nx dx = mesh % x_node(i + 1, j) - mesh % x_node(i, j) dy = mesh % y_node(i + 1, j) - mesh % y_node(i, j) mesh % s_eta(1, i, j) = -dy mesh % s_eta(2, i, j) = dx if (mesh % s_eta(1, i, j)**2 + mesh % s_eta(2, i, j)**2 <= 0.0_wp) then is_ok = .false. write (message, '(A,I0,A,I0,A)') 'mesh_2d: zero-length eta-face at (', i, ',', j, ')' return end if end do end do ! Cell areas (diagonal cross-product) and centres (corner mean). do j = 1, ny do i = 1, nx mesh % vol(i, j) = 0.5_wp * ( & (mesh % x_node(i + 1, j + 1) - mesh % x_node(i, j)) * & (mesh % y_node(i, j + 1) - mesh % y_node(i + 1, j)) - & (mesh % x_node(i, j + 1) - mesh % x_node(i + 1, j)) * & (mesh % y_node(i + 1, j + 1) - mesh % y_node(i, j))) if (mesh % vol(i, j) <= 0.0_wp) then is_ok = .false. write (message, '(A,I0,A,I0,A)') 'mesh_2d: non-positive cell area at (', & i, ',', j, ') — check node ordering / folded grid' return end if mesh % xc(i, j) = 0.25_wp * (mesh % x_node(i, j) + mesh % x_node(i + 1, j) + & mesh % x_node(i + 1, j + 1) + mesh % x_node(i, j + 1)) mesh % yc(i, j) = 0.25_wp * (mesh % y_node(i, j) + mesh % y_node(i + 1, j) + & mesh % y_node(i + 1, j + 1) + mesh % y_node(i, j + 1)) end do end do end subroutine compute_metrics_2d !> Stencil width for the high-order conservative-flux metric operator on n !! nodes: 6 (6th-order) when n>=6, 4 (4th-order) when n>=4, else 2 (the !! classical central scheme). Even widths keep the interior operator a !! symmetric central difference, whose conservative flux telescopes cleanly. pure integer function cons_width(n) result(w) integer, intent(in) :: n if (n >= 6) then w = 6 else if (n >= 4) then w = 4 else w = 2 end if end function cons_width !> Lagrange weights c(1:w) for interpolating/extrapolating a unit-spacing line !! sampled at local positions 0,1,...,w-1 to the (possibly out-of-range) target !! offset t. Used to build the ghost values that close the symmetric stencil at !! domain boundaries. pure subroutine lagrange_unit_weights(w, t, c) integer, intent(in) :: w real(wp), intent(in) :: t real(wp), intent(out) :: c(:) integer :: m, k real(wp) :: p do m = 0, w - 1 p = 1.0_wp do k = 0, w - 1 if (k /= m) p = p * (t - real(k, wp)) / real(m - k, wp) end do c(m + 1) = p end do end subroutine lagrange_unit_weights !> High-order freestream-consistent CONSERVATIVE-FLUX face values of a !! unit-spacing nodal line g(1..n): hf(1..n+1), where face iface lies between !! node iface-1 and node iface. The construction guarantees the telescoping !! identity !! !! hf(i+1) - hf(i) = D[g](i) (the matching high-order central derivative) !! !! EXACTLY (to round-off) at every node, interior AND boundary. Because the !! nodal derivatives AND the transverse face metrics are BOTH built from this !! single operator, and central operators on the two computational axes commute !! discretely, the discrete metric identity Delta_xi(S_xi) + Delta_eta(S_eta)=0 !! holds to round-off on ANY map — the high-order freestream-consistent metric !! of Visbal-Gaitonde (2002) / Nonomura-Iizuka-Fujii (2010) (design spec §5.3). !! !! Interior faces use the symmetric width-w conservative flux (w = cons_width): !! * w=6: [1,-8,37,37,-8,1]/60, whose undivided difference is the 6th-order !! central first derivative; !! * w=4: [-1,7,7,-1]/12 (4th-order central); !! * w=2: the classical [1,1]/2 (2nd-order central — the legacy operator). !! Near a boundary the symmetric stencil is closed by high-order polynomial !! extrapolation of the w nearest interior nodes into the w/2 ghost values it !! needs. Extrapolation keeps the operator a FIXED LINEAR stencil (so the !! cross-axis operators still commute, hence freestream is exact) while keeping !! the closure high order; for the smooth body-fitted maps of interest this !! boundary closure does not degrade the global convergence order. pure subroutine cons_flux_1d_unit(g, hf) real(wp), intent(in) :: g(:) real(wp), intent(out) :: hf(:) integer :: n, w, half, iface, ll, s, m real(wp) :: cf(6), cl(6), acc real(wp) :: gp(1 - 3:size(g) + 3) n = size(g) w = cons_width(n) half = w / 2 select case (w) case (6) cf(1:6) = [1.0_wp, -8.0_wp, 37.0_wp, 37.0_wp, -8.0_wp, 1.0_wp] / 60.0_wp case (4) cf(1:4) = [-1.0_wp, 7.0_wp, 7.0_wp, -1.0_wp] / 12.0_wp case default cf(1:2) = [0.5_wp, 0.5_wp] end select ! Interior nodes verbatim; ghost nodes by high-order Lagrange extrapolation of ! the w nearest boundary nodes, so every face can use the centred stencil. do s = 1, n gp(s) = g(s) end do do m = 1, half ! Left ghost at position 1-m, extrapolated from nodes 1..w (local 0..w-1). call lagrange_unit_weights(w, real(-m, wp), cl) acc = 0.0_wp do s = 1, w acc = acc + cl(s) * g(s) end do gp(1 - m) = acc ! Right ghost at position n+m, extrapolated from nodes n-w+1..n. call lagrange_unit_weights(w, real(m + w - 1, wp), cl) acc = 0.0_wp do s = 1, w acc = acc + cl(s) * g(n - w + s) end do gp(n + m) = acc end do ! Centred width-w conservative flux at every face. Face iface sits between ! node iface-1 and node iface; the stencil spans nodes (iface-1)-(half-1) .. ! (iface-1)+half, i.e. local index ll..ll+w-1. do iface = 1, n + 1 ll = iface - 1 - (half - 1) acc = 0.0_wp do s = 1, w acc = acc + cf(s) * gp(ll + s - 1) end do hf(iface) = acc end do end subroutine cons_flux_1d_unit !> Matching high-order first derivative of a unit-spacing nodal line: the !! undivided difference of the conservative-flux face values !! (cons_flux_1d_unit), so dg(i) reproduces the high-order central derivative !! AND stays bit-consistent with the face metrics built from the same operator !! (the property that makes freestream telescope to round-off, design spec §5.3). pure subroutine deriv_1d_unit(g, dg) real(wp), intent(in) :: g(:) real(wp), intent(out) :: dg(:) real(wp) :: hf(size(g) + 1) integer :: i, n n = size(g) call cons_flux_1d_unit(g, hf) do i = 1, n dg(i) = hf(i + 1) - hf(i) end do end subroutine deriv_1d_unit !> Compute NODAL transformation metrics for curvilinear FDM from the node !! coordinates already in mesh%x_node / mesh%y_node (shape (nx, ny), the Plot3D !! points treated as nodes; unit computational spacing dxi=deta=1). !! !! Allocates (with stat=) and fills: !! * jac(nx, ny) = x_xi*y_eta - x_eta*y_xi at NODES (must be > 0) !! * sx_xi(2, nx+1, ny) = (y_eta, -x_eta) at XI-faces !! * sx_eta(2, nx, ny+1) = (-y_xi, x_xi) at ETA-faces !! and sets mesh%fdm_curvilinear = .true. !! !! Geometry derivatives AND the transverse face metrics are built from the !! SAME high-order conservative-flux operator (cons_flux_1d_unit: a symmetric !! width-w central conservative flux — 6th order for n>=6 — closed by high-order !! extrapolation at boundaries). The nodal derivatives are the undivided !! difference of that operator's face values (deriv_1d_unit), and the face !! metric at face iface IS the operator's conservative-flux value of the !! transverse nodal derivative, so the telescoped first difference of the face !! metric reproduces the high-order nodal derivative EXACTLY at EVERY node !! (interior and boundary). Because deriv_xi and deriv_eta are tensor-product !! operators on separate axes they commute discretely, so with this consistent !! HIGH-ORDER evaluation the discrete identities !! Delta_xi(y_eta) - Delta_eta(y_xi) = 0, Delta_xi(x_eta) - Delta_eta(x_xi) = 0 !! hold to machine precision at every node even on a fully non-separable map — !! the constant flow telescopes to round-off (the freestream/GCL gate) WITHOUT !! capping the curved-grid order at 2 (the high-order metric of Visbal-Gaitonde !! 2002 / Nonomura-Iizuka-Fujii 2010, design spec §5.3). !! !! Face index convention: the xi-face between node i and node i+1 is stored at !! index i+1 in 1..nx+1, so index 1 is the left boundary face (node 1) and !! index nx+1 is the right boundary face (node nx). Likewise for eta-faces. subroutine compute_metrics_fdm_2d(mesh, is_ok, message) type(mesh_2d_t), intent(inout) :: mesh logical, intent(out) :: is_ok character(len=*), intent(out) :: message integer :: nx, ny, i, j, astat real(wp), allocatable :: x_xi(:, :), x_eta(:, :), y_xi(:, :), y_eta(:, :) real(wp), allocatable :: line(:), dline(:), hfx(:), hfy(:) is_ok = .true. message = '' nx = mesh % nx; ny = mesh % ny if (nx < 3 .or. ny < 3) then is_ok = .false. message = 'mesh_2d: curvilinear FDM requires nx >= 3 and ny >= 3 nodes per axis' return end if if (.not. (allocated(mesh % x_node) .and. allocated(mesh % y_node))) then is_ok = .false. message = 'mesh_2d: compute_metrics_fdm_2d requires x_node/y_node' return end if allocate (x_xi(nx, ny), x_eta(nx, ny), y_xi(nx, ny), y_eta(nx, ny), stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: allocation failed for nodal derivative scratch'; return end if ! xi-derivatives (along i) row by row. allocate (line(nx), dline(nx), stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: allocation failed for xi-derivative line'; return end if do j = 1, ny line = mesh % x_node(:, j); call deriv_1d_unit(line, dline); x_xi(:, j) = dline line = mesh % y_node(:, j); call deriv_1d_unit(line, dline); y_xi(:, j) = dline end do deallocate (line, dline, stat=astat) if (astat /= 0) error stop 'mesh_2d: xi-derivative line deallocate failed' ! eta-derivatives (along j) column by column. allocate (line(ny), dline(ny), stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: allocation failed for eta-derivative line'; return end if do i = 1, nx line = mesh % x_node(i, :); call deriv_1d_unit(line, dline); x_eta(i, :) = dline line = mesh % y_node(i, :); call deriv_1d_unit(line, dline); y_eta(i, :) = dline end do deallocate (line, dline, stat=astat) if (astat /= 0) error stop 'mesh_2d: eta-derivative line deallocate failed' ! Allocate the metric arrays (re-allocate defensively on a dimension change). ! Deallocate any stale array with checked stat= (do not silently overwrite the ! shared stat between the defensive deallocates and the re-allocate). if (allocated(mesh % jac)) then deallocate (mesh % jac, stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: deallocation of stale jac failed'; return end if end if if (allocated(mesh % sx_xi)) then deallocate (mesh % sx_xi, stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: deallocation of stale sx_xi failed'; return end if end if if (allocated(mesh % sx_eta)) then deallocate (mesh % sx_eta, stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: deallocation of stale sx_eta failed'; return end if end if allocate (mesh % jac(nx, ny), mesh % sx_xi(2, nx + 1, ny), mesh % sx_eta(2, nx, ny + 1), stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: allocation failed for nodal metric arrays'; return end if ! Jacobian at nodes (guard J > 0). do j = 1, ny do i = 1, nx mesh % jac(i, j) = x_xi(i, j) * y_eta(i, j) - x_eta(i, j) * y_xi(i, j) if (mesh % jac(i, j) <= 0.0_wp) then is_ok = .false. write (message, '(A,I0,A,I0,A)') 'mesh_2d: non-positive Jacobian at node (', & i, ',', j, ') — check node ordering / folded grid' return end if end do end do allocate (hfx(nx + 1), hfy(ny + 1), stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: allocation failed for conservative-flux face scratch'; return end if ! xi-face metrics S_xi = (y_eta, -x_eta): the high-order conservative-flux ! interpolation (cons_flux_1d_unit) of the transverse nodal derivative ALONG i, ! evaluated at faces 1..nx+1. Its telescoped difference is the high-order ! nodal xi-derivative, so freestream telescopes to round-off (see header). do j = 1, ny call cons_flux_1d_unit(y_eta(:, j), hfx); mesh % sx_xi(1, :, j) = hfx call cons_flux_1d_unit(x_eta(:, j), hfx); mesh % sx_xi(2, :, j) = -hfx end do ! eta-face metrics S_eta = (-y_xi, x_xi): the same operator ALONG j, evaluated ! at faces 1..ny+1. do i = 1, nx call cons_flux_1d_unit(y_xi(i, :), hfy); mesh % sx_eta(1, i, :) = -hfy call cons_flux_1d_unit(x_xi(i, :), hfy); mesh % sx_eta(2, i, :) = hfy end do mesh % fdm_curvilinear = .true. deallocate (x_xi, x_eta, y_xi, y_eta, hfx, hfy, stat=astat) if (astat /= 0) error stop 'mesh_2d: nodal derivative scratch deallocate failed' end subroutine compute_metrics_fdm_2d !> Build a curvilinear-FDM mesh from a single-block Plot3D grid: the Plot3D !! points are treated as NODES (mesh % nx = ni, mesh % ny = nj, nodes lying ON !! the boundaries), mesh % uniform = .false., dx/dy = 0 (unused on the !! curvilinear path), then compute_metrics_fdm_2d fills the nodal metrics. !! Distinct from build_mesh_2d_global's 'plot3d' branch, which treats the points !! as CELL CORNERS (nx = ni-1 cells) and builds finite-volume metrics. subroutine build_mesh_2d_fdm_curvilinear(mesh, grid_file, grid_binary, is_ok, message) type(mesh_2d_t), intent(inout) :: mesh character(len=*), intent(in) :: grid_file logical, intent(in) :: grid_binary logical, intent(out) :: is_ok character(len=*), intent(out) :: message integer :: ni, nj, astat real(wp), allocatable :: xg(:, :), yg(:, :) call read_plot3d_2d(grid_file, grid_binary, xg, yg, ni, nj, is_ok, message) if (.not. is_ok) return mesh % uniform = .false. mesh % nx = ni mesh % ny = nj mesh % dx_uniform = 0.0_wp mesh % dy_uniform = 0.0_wp call move_alloc(xg, mesh % x_node) call move_alloc(yg, mesh % y_node) ! Drop any stale finite-volume metric arrays so a re-used mesh_2d_t cannot ! carry FVM metrics into the curvilinear-FDM residual. if (allocated(mesh % vol)) deallocate (mesh % vol, stat=astat) if (allocated(mesh % s_xi)) deallocate (mesh % s_xi, stat=astat) if (allocated(mesh % s_eta)) deallocate (mesh % s_eta, stat=astat) if (allocated(mesh % xc)) deallocate (mesh % xc, stat=astat) if (allocated(mesh % yc)) deallocate (mesh % yc, stat=astat) call compute_metrics_fdm_2d(mesh, is_ok, message) end subroutine build_mesh_2d_fdm_curvilinear !> Build a uniform Cartesian NODAL mesh for the finite-difference (FDM) path. !! Unlike build_mesh_2d_global's 'uniform' branch (which stores CELL counts and !! derives FV metrics), this stores NODE counts directly: mesh % nx = npx, !! mesh % ny = npy, with grid points lying ON both boundaries. The spacing is !! dx = (x_right - x_left) / (npx - 1) (equal to the FV cell size for the same !! number of cells). The FV metric arrays (vol, s_xi, s_eta, xc, yc) are left !! UNALLOCATED on purpose — the uniform FDM residual divides by the scalar !! dx/dy, not by vol — so this builder does NOT call compute_metrics_2d. subroutine build_mesh_2d_uniform_nodal(mesh, npx, npy, x_left, x_right, y_left, y_right) type(mesh_2d_t), intent(inout) :: mesh integer, intent(in) :: npx, npy real(wp), intent(in) :: x_left, x_right, y_left, y_right integer :: i, j, astat if (npx < 2 .or. npy < 2) & error stop 'mesh_2d: build_mesh_2d_uniform_nodal requires npx >= 2 and npy >= 2' mesh % uniform = .true. mesh % nx = npx mesh % ny = npy mesh % dx_uniform = (x_right - x_left) / real(npx - 1, wp) mesh % dy_uniform = (y_right - y_left) / real(npy - 1, wp) if (allocated(mesh % x_node)) then deallocate (mesh % x_node, stat=astat) if (astat /= 0) error stop 'mesh_2d: x_node deallocate failed' end if if (allocated(mesh % y_node)) then deallocate (mesh % y_node, stat=astat) if (astat /= 0) error stop 'mesh_2d: y_node deallocate failed' end if ! FV metrics are left UNALLOCATED on purpose for the nodal mesh. If this ! mesh_2d_t was previously initialised via build_mesh_2d_global + ! compute_metrics_2d, stale vol/s_xi/s_eta/xc/yc would remain allocated with ! the wrong dimensions. Deallocate them so the nodal mesh is truly metric-free. if (allocated(mesh % vol)) then deallocate (mesh % vol, stat=astat) if (astat /= 0) error stop 'mesh_2d: vol deallocate failed' end if if (allocated(mesh % s_xi)) then deallocate (mesh % s_xi, stat=astat) if (astat /= 0) error stop 'mesh_2d: s_xi deallocate failed' end if if (allocated(mesh % s_eta)) then deallocate (mesh % s_eta, stat=astat) if (astat /= 0) error stop 'mesh_2d: s_eta deallocate failed' end if if (allocated(mesh % xc)) then deallocate (mesh % xc, stat=astat) if (astat /= 0) error stop 'mesh_2d: xc deallocate failed' end if if (allocated(mesh % yc)) then deallocate (mesh % yc, stat=astat) if (astat /= 0) error stop 'mesh_2d: yc deallocate failed' end if allocate (mesh % x_node(npx, npy), mesh % y_node(npx, npy), stat=astat) if (astat /= 0) error stop 'mesh_2d: allocation failed for nodal node arrays' do j = 1, npy do i = 1, npx mesh % x_node(i, j) = x_left + mesh % dx_uniform * real(i - 1, wp) mesh % y_node(i, j) = y_left + mesh % dy_uniform * real(j - 1, wp) end do end do end subroutine build_mesh_2d_uniform_nodal !> Slice the (global) metric arrays of `mesh` into this rank's local sub-block. !! Cells ix_first..ix_first+nx_local-1, iy_first..iy_first+ny_local-1 (1-based !! global). xi-faces span nx_local+1 columns starting at global column !! ix_first; eta-faces span ny_local+1 rows starting at iy_first. Shared !! inter-rank faces are identical on both ranks because both slice the same !! global metrics. No-op for uniform grids (they use scalar dx/dy). subroutine build_mesh_2d_local(mesh, ix_first, iy_first, nx_local, ny_local, is_ok, message) type(mesh_2d_t), intent(inout) :: mesh integer, intent(in) :: ix_first, iy_first, nx_local, ny_local logical, intent(out) :: is_ok character(len=*), intent(out) :: message real(wp), allocatable :: vol_l(:, :), xc_l(:, :), yc_l(:, :) real(wp), allocatable :: sxi_l(:, :, :), seta_l(:, :, :) real(wp), allocatable :: xn_l(:, :), yn_l(:, :), jac_l(:, :) integer :: i, j, astat, i0, j0 is_ok = .true. message = '' if (mesh % uniform) return ! uniform path uses scalar dx/dy; metrics unused if (ix_first < 1 .or. iy_first < 1 .or. & ix_first + nx_local - 1 > mesh % nx .or. iy_first + ny_local - 1 > mesh % ny) then is_ok = .false. message = 'mesh_2d: build_mesh_2d_local index range out of global bounds' return end if i0 = ix_first - 1; j0 = iy_first - 1 if (mesh % fdm_curvilinear) then ! Curvilinear-FDM nodal-metric slice. The Plot3D points are NODES, so ! mesh % nx/ny are NODE counts and the local block owns nx_local x ny_local ! NODES at global (ix_first, iy_first). xi-faces number nx+1 (sx_xi), so the ! local block has nx_local+1 xi-faces starting at global face ix_first; ! eta-faces number ny+1. Both ranks slice the SAME global metric, so a ! shared inter-rank face is bit-identical on the two ranks (the FVM ! curvilinear MPI invariant, applied to the nodal metrics). Distinct from ! the FVM branch below, whose x_node array is nx+1 (cell corners); here ! x_node/y_node are node-count arrays sliced to (nx_local, ny_local). allocate (jac_l(nx_local, ny_local), & sxi_l(2, nx_local + 1, ny_local), seta_l(2, nx_local, ny_local + 1), & xn_l(nx_local, ny_local), yn_l(nx_local, ny_local), stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: build_mesh_2d_local (fdm) allocation failed'; return end if do j = 1, ny_local do i = 1, nx_local jac_l(i, j) = mesh % jac(i0 + i, j0 + j) xn_l(i, j) = mesh % x_node(i0 + i, j0 + j) yn_l(i, j) = mesh % y_node(i0 + i, j0 + j) end do end do do j = 1, ny_local do i = 1, nx_local + 1 sxi_l(:, i, j) = mesh % sx_xi(:, i0 + i, j0 + j) end do end do do j = 1, ny_local + 1 do i = 1, nx_local seta_l(:, i, j) = mesh % sx_eta(:, i0 + i, j0 + j) end do end do call move_alloc(jac_l, mesh % jac) call move_alloc(sxi_l, mesh % sx_xi) call move_alloc(seta_l, mesh % sx_eta) call move_alloc(xn_l, mesh % x_node) call move_alloc(yn_l, mesh % y_node) mesh % nx = nx_local mesh % ny = ny_local return end if allocate (vol_l(nx_local, ny_local), xc_l(nx_local, ny_local), yc_l(nx_local, ny_local), & sxi_l(2, nx_local + 1, ny_local), seta_l(2, nx_local, ny_local + 1), & xn_l(nx_local + 1, ny_local + 1), yn_l(nx_local + 1, ny_local + 1), stat=astat) if (astat /= 0) then is_ok = .false.; message = 'mesh_2d: build_mesh_2d_local allocation failed'; return end if do j = 1, ny_local do i = 1, nx_local vol_l(i, j) = mesh % vol(i0 + i, j0 + j) xc_l(i, j) = mesh % xc(i0 + i, j0 + j) yc_l(i, j) = mesh % yc(i0 + i, j0 + j) end do end do ! Local node block: nodes ix_first..ix_first+nx_local (nx_local+1 of them). do j = 1, ny_local + 1 do i = 1, nx_local + 1 xn_l(i, j) = mesh % x_node(i0 + i, j0 + j) yn_l(i, j) = mesh % y_node(i0 + i, j0 + j) end do end do do j = 1, ny_local do i = 1, nx_local + 1 sxi_l(:, i, j) = mesh % s_xi(:, i0 + i, j0 + j) end do end do do j = 1, ny_local + 1 do i = 1, nx_local seta_l(:, i, j) = mesh % s_eta(:, i0 + i, j0 + j) end do end do call move_alloc(vol_l, mesh % vol) call move_alloc(xc_l, mesh % xc) call move_alloc(yc_l, mesh % yc) call move_alloc(sxi_l, mesh % s_xi) call move_alloc(seta_l, mesh % s_eta) call move_alloc(xn_l, mesh % x_node) call move_alloc(yn_l, mesh % y_node) mesh % nx = nx_local mesh % ny = ny_local end subroutine build_mesh_2d_local end module mesh_2d