diff --git a/.github/workflows/fprettify-lint.yml b/.github/workflows/fprettify-lint.yml index b002184a6..8810bf167 100644 --- a/.github/workflows/fprettify-lint.yml +++ b/.github/workflows/fprettify-lint.yml @@ -5,11 +5,13 @@ on: branches: - main paths-ignore: + - 'docs/' - 'AUTHORS.md' - 'LICENSE.md' - 'README.md' pull_request: paths-ignore: + - 'docs/' - 'AUTHORS.md' - 'LICENSE.md' - 'README.md' diff --git a/.github/workflows/main-docs.yml b/.github/workflows/main-docs.yml index 8fd5c718c..16f4387d4 100644 --- a/.github/workflows/main-docs.yml +++ b/.github/workflows/main-docs.yml @@ -3,6 +3,10 @@ on: push: branches: - main + paths: + - 'docs/**' + - 'mkdocs.yml' + - 'self.md' jobs: build: diff --git a/docs/Learning/ProvableStability.md b/docs/Learning/ProvableStability.md index 1bcc58396..a15b8455d 100644 --- a/docs/Learning/ProvableStability.md +++ b/docs/Learning/ProvableStability.md @@ -98,7 +98,28 @@ where $s$ is the solution and $u$ is a velocity field. A convex entropy function ## Split form equations +Split-form equations rewrite the flux divergence so that the resulting discretization satisfies a discrete entropy conservation law. For a scalar conservation law $s_t + (u s)_x = 0$, we can use the product rule to write the equivalent split form: +\begin{equation} +s_t + \frac{1}{2}\left( (u s)_x + u s_x + s u_x \right) = 0 +\end{equation} + +Discretizing the split form with the SBP derivative matrix $D$ and a symmetric two-point flux $f^\#(s_i, s_n)$ leads to the split-form DGSEM: + +\begin{equation} +\frac{ds_i}{dt} = -2 \sum_n D_{\text{split},n,i} \, f^\#(s_i, s_n) - \frac{1}{w_i} \left( f_{\text{Riemann}} - f_{\text{local}} \right)\bigg|_{\partial \Omega_e} +\end{equation} + +where $D_\text{split} = D - \frac{1}{2}M^{-1}B$ is the skew-symmetric split-form derivative operator (see [Split-Form DGSEM](SplitFormDGSEM.md) for details). The skew-symmetry of $D_\text{split}$ ensures the volume integral contributes zero to the discrete entropy rate; all entropy change passes through the surface Riemann flux. ## Two-point flux +A two-point flux $\mathbf{f}^\#(\mathbf{s}_L, \mathbf{s}_R)$ is a symmetric, consistent numerical flux that satisfies the Tadmor entropy conservation condition: + +\begin{equation} +(\mathbf{w}_R - \mathbf{w}_L)^T \mathbf{f}^\# = \Psi_R - \Psi_L +\end{equation} + +where $\mathbf{w} = \partial \eta / \partial \mathbf{s}$ are the entropy variables and $\Psi$ is the entropy flux potential. When used in the split-form volume integral, this condition guarantees that the discrete volume contribution to $d\eta/dt$ vanishes identically. + +For linear advection with $f = as$ and $\eta = s^2/2$, the arithmetic mean $f^\#(s_L, s_R) = a(s_L + s_R)/2$ satisfies this condition. For nonlinear systems (Euler, shallow water), deriving the entropy-conserving flux is more involved; see the references in [Split-Form DGSEM](SplitFormDGSEM.md). diff --git a/docs/Learning/SoftwareArchitecture.md b/docs/Learning/SoftwareArchitecture.md index a675a27eb..c81e2e5bc 100644 --- a/docs/Learning/SoftwareArchitecture.md +++ b/docs/Learning/SoftwareArchitecture.md @@ -78,7 +78,7 @@ and in 3-D F(0:N,0:N,0:N,0:N,1:nEl,1:nvar,1:3) ``` -The first dimension refers to the directional component (physical or computational) of the vector. The second dimension is a computational coordinate direction and variations in this dimensions are equivalent to looping over $n$ in the two-point vector representation shown above. The remaining array dimensions of size `0:N` are computational coordinate directions and are the usual `(i,j,k)` dimensions, in this order. This is followed by the dimension over the solution variables and then the elements. +The first dimension is a computational coordinate direction and variations in this dimension are equivalent to looping over $n$ in the two-point vector representation shown above. The next dimensions are `(i,j)` (2-D) or `(i,j,k)` (3-D) computational coordinate directions, followed by the elements, solution variables, and the directional component of the vector `(...,nEl,nvar,idir)`, respectively. ### Tensors diff --git a/docs/Learning/SplitFormDGSEM.md b/docs/Learning/SplitFormDGSEM.md new file mode 100644 index 000000000..f26f72e76 --- /dev/null +++ b/docs/Learning/SplitFormDGSEM.md @@ -0,0 +1,166 @@ +# Split-Form Entropy-Conserving DGSEM + +This page documents the split-form discontinuous Galerkin spectral element method (DGSEM) implemented in SELF's `ECDGModel2D_t` and `ECDGModel3D_t` classes. The formulation follows Gassner, Winters, and Kopriva (2016) and is aligned with the [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) implementation for curved meshes. + +## Mathematical Background + +### Conservation Law +We consider a system of conservation laws on a domain $\Omega$: + +$$ +\mathbf{s}_t + \vec{\nabla} \cdot \vec{\mathbf{f}}(\mathbf{s}) = \mathbf{q} +$$ + +where $\mathbf{s}$ is the state vector, $\vec{\mathbf{f}}$ is the flux, and $\mathbf{q}$ is a source term. + +### Entropy Conservation +A convex entropy function $\eta(\mathbf{s})$ satisfies the entropy conservation law $\eta_t + \vec{\nabla} \cdot \vec{\Psi} = 0$ (see [Provable Stability](ProvableStability.md)). A numerical scheme is **entropy-conserving** if it satisfies a discrete analogue of this conservation law, and **entropy-stable** if entropy can only decrease (mimicking physical dissipation at discontinuities). + +### Two-Point Entropy-Conserving Flux +An entropy-conserving two-point flux $\mathbf{f}^\#(\mathbf{s}_L, \mathbf{s}_R)$ is a symmetric function of two states satisfying the Tadmor condition: + +$$ +(\mathbf{w}_R - \mathbf{w}_L)^T \mathbf{f}^\#(\mathbf{s}_L, \mathbf{s}_R) = \Psi_R - \Psi_L +$$ + +where $\mathbf{w} = \partial \eta / \partial \mathbf{s}$ are the entropy variables and $\Psi$ is the entropy flux potential. + +**Example:** For linear advection $f = a s$ with entropy $\eta = s^2/2$, the arithmetic mean $f^\#(s_L, s_R) = a(s_L + s_R)/2$ is entropy-conserving. + +## The Split-Form Derivative Operator + +### Standard Derivative Matrix $D$ +The standard DGSEM derivative matrix $D$ satisfies $D \cdot \mathbf{1} = 0$ (derivative of a constant is zero) and the SBP property: + +$$ +M D + D^T M = B +$$ + +where $M = \text{diag}(w_0, \ldots, w_N)$ is the mass matrix (quadrature weights) and $B$ is the boundary operator. + +### Split-Form Matrix $D_\text{split}$ +The split-form derivative matrix is defined as: + +$$ +D_\text{split} = D - \tfrac{1}{2} M^{-1} B +$$ + +**Key property:** $D_\text{split}$ is skew-symmetric under the $M$ inner product: + +$$ +M D_\text{split} + D_\text{split}^T M = 0 +$$ + +This follows directly from the SBP property: $M D_\text{split} + D_\text{split}^T M = (M D + D^T M) - B = B - B = 0$. + +**$D_\text{split}$ is NOT an SBP operator.** Its weighted symmetric part is zero (not the boundary term $B$). This means the volume integral formed with $D_\text{split}$ contributes **zero** to the entropy rate. All entropy change passes through the surface term. + +For Gauss-Lobatto-Legendre (GLL) nodes, $D_\text{split}$ reduces to $(D - D^T)/2$, and SELF requires GLL nodes for `TwoPointVector` types. + +In SELF, the split-form matrix is stored as `interp%dSplitMatrix` in the Lagrange class: + +```fortran +dSplitMatrix(ii,i) = dMatrix(ii,i) & + - 0.5*(bMatrix(i,2)*bMatrix(ii,2) - bMatrix(i,1)*bMatrix(ii,1)) / qWeights(i) +``` + +### Relationship to Trixi.jl +Trixi.jl defines `derivative_split = 2D - M^{-1}B` and applies it **without** a factor of 2. SELF defines `dSplitMatrix = D - 0.5 M^{-1}B` and multiplies by 2 in the divergence formula. These are algebraically identical: + +$$ +2 \cdot D_\text{split}^\text{SELF} = D_\text{split}^\text{Trixi} +$$ + + +## The EC-DGSEM Semidiscretization + +### Full Tendency +The `ECDGModel2D_t` computes the time tendency as: + +$$ +\frac{d\mathbf{s}}{dt} = \mathbf{q} - \frac{1}{J} \left[ 2 \sum_n D_{\text{split},n,i} \, \tilde{F}^r(n,i,j) \;+\; M^{-1} B^T \mathbf{f}_\text{Riemann} \right] +$$ + +where: + +- $\tilde{F}^r(n,i,j)$ is a **scalar contravariant** two-point flux for direction $r$ +- $\mathbf{f}_\text{Riemann}$ is the surface Riemann flux (e.g., Local Lax-Friedrichs) +- $J$ is the element Jacobian + +### Equivalence to the Strong-Form Penalty +Using $D_\text{split}$ for the volume combined with the plain Riemann flux on the surface is algebraically identical to using the standard $D$ for the volume with a **penalty** $(f_\text{Riemann} - f_\text{local})$ on the surface: + +$$ +2 D_\text{split} \tilde{F} + M^{-1} B^T f_\text{Riemann} = 2 D \tilde{F} + M^{-1} B^T (f_\text{Riemann} - \tilde{F}_\text{boundary}) +$$ + +SELF uses the left-hand form (matching Trixi.jl's `SurfaceIntegralWeakForm`). + +## Contravariant Two-Point Fluxes on Curved Meshes + +### Convention +Following Trixi.jl for curved meshes, `interior(n,i,j,iEl,iVar,r)` stores a **scalar** contravariant two-point flux for the $r$-th computational direction. Each direction uses its own node pairing and its own averaged metric: + +| Direction $r$ | Node pair | Averaged metric | +|:---:|---|---| +| $\xi^1$ | $(i,j) \leftrightarrow (n,j)$ | $\overline{J\mathbf{a}^1} = \tfrac{1}{2}\left(J\mathbf{a}^1_{i,j} + J\mathbf{a}^1_{n,j}\right)$ | +| $\xi^2$ | $(i,j) \leftrightarrow (i,n)$ | $\overline{J\mathbf{a}^2} = \tfrac{1}{2}\left(J\mathbf{a}^2_{i,j} + J\mathbf{a}^2_{i,n}\right)$ | +| $\xi^3$ | $(i,j,k) \leftrightarrow (i,j,n)$ | $\overline{J\mathbf{a}^3} = \tfrac{1}{2}\left(J\mathbf{a}^3_{i,j,k} + J\mathbf{a}^3_{i,j,n}\right)$ | + +The contravariant flux is the dot product of the averaged metric with the physical EC flux: + +$$ +\tilde{F}^r(n,i,j) = \sum_d \overline{J a^r_d} \; f^\#_d(\mathbf{s}_L, \mathbf{s}_R) +$$ + +where $f^\#_d$ is the $d$-th physical component of the entropy-conserving two-point flux evaluated between the direction-specific node pair. + +### Why Not Physical-Space Storage? +A single `interior(n,i,j,...,d)` array with $d$ indexing physical directions cannot correctly serve both the $\xi^1$ and $\xi^2$ sums, because the same index $n$ refers to different physical node partners in each direction. Storing pre-projected scalar contravariant fluxes (one per computational direction) avoids this cross-contamination entirely. + +### Implementation +`TwoPointFluxMethod` in `ECDGModel2D_t` computes the contravariant projection: + +```fortran +! xi^1: pair (i,j)-(nn,j), project onto avg(Ja^1) +sR = solution(nn,j,iel,:) +Fphys = twopointflux2d(sL, sR) ! physical EC flux (nvar x 2) +Fc = sum_d avg(Ja^1_d) * Fphys(:,d) ! scalar contravariant +twoPointFlux%interior(nn,i,j,...,1) = Fc + +! xi^2: pair (i,j)-(i,nn), project onto avg(Ja^2) +sR = solution(i,nn,iel,:) +Fphys = twopointflux2d(sL, sR) +Fc = sum_d avg(Ja^2_d) * Fphys(:,d) +twoPointFlux%interior(nn,i,j,...,2) = Fc +``` + +`MappedDivergence` then simply applies `Divergence / J` -- no metric operations inside. + +## Implementing a Concrete EC Model + +To create an entropy-conserving model, extend `ECDGModel2D_t` (or `ECDGModel3D_t`) and override: + +| Procedure | Purpose | +|-----------|---------| +| `SetNumberOfVariables` | Set `this%nvar` | +| `twopointflux2d(sL, sR)` | Return the physical EC two-point flux (nvar x 2 array) | +| `riemannflux2d(sL, sR, dsdx, nhat)` | Return the surface Riemann flux (nvar array). Use LLF for symmetric dissipation. | +| `entropy_func(s)` | Return the scalar entropy $\eta(\mathbf{s})$ for diagnostics | +| `SetMetadata` | (Optional) name and units for each variable | + +The base class handles the split-form volume integral, metric averaging, surface correction, time integration (inherited from `DGModel2D_t`), and I/O. + +**Example:** `ECAdvection2D_t` implements entropy-conserving linear advection with: +- `twopointflux2d`: arithmetic mean $f^\# = \mathbf{a}(s_L + s_R)/2$ +- `riemannflux2d`: Local Lax-Friedrichs $f_\text{LLF} = \tfrac{1}{2}\left[u_n(s_L+s_R) - |\mathbf{a}|(s_R-s_L)\right]$ +- `entropy_func`: $\eta = s^2/2$ + + +## References + +1. Gassner, G. J., Winters, A. R., & Kopriva, D. A. (2016). Split form nodal discontinuous Galerkin schemes with summation-by-parts property for the compressible Euler equations. *Journal of Computational Physics*, 327, 39-66. + +2. Winters, A. R., Kopriva, D. A., Gassner, G. J., & Hindenlang, F. (2021). Construction of Modern Robust Nodal Discontinuous Galerkin Spectral Element Methods for the Compressible Navier-Stokes Equations. *Lecture Notes in Computational Science and Engineering*, Springer. + +3. Ranocha, H., Schlottke-Lakemper, M., Winters, A. R., Faulhaber, E., Chan, J., & Gassner, G. J. (2022). Adaptive numerical simulations with Trixi.jl: A case study of Julia for scientific computing. *Proceedings of the JuliaCon Conferences*. diff --git a/src/SELF_DGModel1D_t.f90 b/src/SELF_DGModel1D_t.f90 index 26256bb21..75122257d 100644 --- a/src/SELF_DGModel1D_t.f90 +++ b/src/SELF_DGModel1D_t.f90 @@ -307,7 +307,7 @@ subroutine CalculateEntropy_DGModel1D_t(this) do i = 1,this%solution%interp%N+1 J = this%geometry%dxds%interior(i,iel,1) s(1:this%solution%nvar) = this%solution%interior(i,iel,1:this%solution%nvar) - e = e+this%entropy_func(s)*J + e = e+this%entropy_func(s)*J*this%solution%interp%qWeights(i) enddo enddo diff --git a/src/SELF_DGModel2D_t.f90 b/src/SELF_DGModel2D_t.f90 index df9751af6..d87b19a05 100644 --- a/src/SELF_DGModel2D_t.f90 +++ b/src/SELF_DGModel2D_t.f90 @@ -343,7 +343,9 @@ subroutine CalculateEntropy_DGModel2D_t(this) do i = 1,this%solution%interp%N+1 jac = abs(this%geometry%J%interior(i,j,iel,1)) s = this%solution%interior(i,j,iel,1:this%nvar) - e = e+this%entropy_func(s)*jac + e = e+this%entropy_func(s)*jac* & + this%solution%interp%qWeights(i)* & + this%solution%interp%qWeights(j) enddo enddo enddo diff --git a/src/SELF_DGModel3D_t.f90 b/src/SELF_DGModel3D_t.f90 index 6b8651163..f0d9c3de1 100644 --- a/src/SELF_DGModel3D_t.f90 +++ b/src/SELF_DGModel3D_t.f90 @@ -344,7 +344,10 @@ subroutine CalculateEntropy_DGModel3D_t(this) do i = 1,this%solution%interp%N+1 jac = abs(this%geometry%J%interior(i,j,k,iel,1)) s = this%solution%interior(i,j,k,iel,1:this%nvar) - e = e+this%entropy_func(s)*jac + e = e+this%entropy_func(s)*jac* & + this%solution%interp%qWeights(i)* & + this%solution%interp%qWeights(j)* & + this%solution%interp%qWeights(k) enddo enddo enddo diff --git a/src/SELF_ECAdvection2D_t.f90 b/src/SELF_ECAdvection2D_t.f90 new file mode 100644 index 000000000..807df153f --- /dev/null +++ b/src/SELF_ECAdvection2D_t.f90 @@ -0,0 +1,131 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECAdvection2D_t + !! Entropy-conserving linear scalar advection equation in 2-D. + !! + !! Solves: + !! du/dt + u_a * du/dx + v_a * du/dy = 0 + !! + !! Entropy function: eta(u) = u^2 / 2 (same as Trixi.jl) + !! + !! EC two-point flux (arithmetic mean, entropy-conserving for eta = u^2/2): + !! F^EC(uL, uR) = (u_a * (uL+uR)/2, v_a * (uL+uR)/2) + !! + !! Surface Riemann flux (Local Lax-Friedrichs / Rusanov): + !! F_Riemann = 0.5 * (a.n) * (uL+uR) - 0.5 * |a| * (uR-uL) + !! where |a| = sqrt(u^2 + v^2) is the maximum wave speed. + !! This is dissipative (entropy-stable) and reduces to the central flux + !! when uL = uR (no dissipation at no-normal-flow or mirror boundaries). + + use SELF_ECDGModel2D_t + use SELF_mesh + + implicit none + + type,extends(ECDGModel2D_t) :: ECAdvection2D_t + + real(prec) :: u ! x-component of advection velocity + real(prec) :: v ! y-component of advection velocity + + contains + + procedure :: entropy_func => entropy_func_ECAdvection2D_t + procedure :: twopointflux2d => twopointflux2d_ECAdvection2D_t + procedure :: riemannflux2d => riemannflux2d_ECAdvection2D_t + + ! Mirror BC: sets extBoundary = interior boundary state so that + ! sR = sL at all domain faces. Used to eliminate boundary dissipation + ! when testing entropy conservation. + procedure :: hbc2d_NoNormalFlow => hbc2d_NoNormalFlow_ECAdvection2D_t + + endtype ECAdvection2D_t + +contains + + pure function entropy_func_ECAdvection2D_t(this,s) result(e) + !! Quadratic entropy: eta(u) = u^2 / 2 + class(ECAdvection2D_t),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec) :: e + + e = 0.5_prec*s(1)*s(1) + if(.false.) e = e+this%u ! suppress unused-dummy-argument warning + + endfunction entropy_func_ECAdvection2D_t + + pure function twopointflux2d_ECAdvection2D_t(this,sL,sR) result(flux) + !! Arithmetic-mean two-point flux for linear advection. + !! Entropy-conserving with respect to eta(u) = u^2/2. + class(ECAdvection2D_t),intent(in) :: this + real(prec),intent(in) :: sL(1:this%nvar) + real(prec),intent(in) :: sR(1:this%nvar) + real(prec) :: flux(1:this%nvar,1:2) + ! Local + real(prec) :: savg + + savg = 0.5_prec*(sL(1)+sR(1)) + flux(1,1) = this%u*savg + flux(1,2) = this%v*savg + + endfunction twopointflux2d_ECAdvection2D_t + + pure function riemannflux2d_ECAdvection2D_t(this,sL,sR,dsdx,nhat) result(flux) + !! Local Lax-Friedrichs (Rusanov) Riemann flux for linear advection. + !! Entropy-stable: provides symmetric dissipation at element interfaces. + !! lambda_max = sqrt(u^2 + v^2) is the maximum wave speed. + class(ECAdvection2D_t),intent(in) :: this + real(prec),intent(in) :: sL(1:this%nvar) + real(prec),intent(in) :: sR(1:this%nvar) + real(prec),intent(in) :: dsdx(1:this%nvar,1:2) + real(prec),intent(in) :: nhat(1:2) + real(prec) :: flux(1:this%nvar) + ! Local + real(prec) :: un,lam + + un = this%u*nhat(1)+this%v*nhat(2) + lam = sqrt(this%u*this%u+this%v*this%v) + flux(1) = 0.5_prec*(un*(sL(1)+sR(1))-lam*(sR(1)-sL(1))) + if(.false.) flux(1) = flux(1)+dsdx(1,1) ! suppress unused-dummy-argument warning + + endfunction riemannflux2d_ECAdvection2D_t + + pure function hbc2d_NoNormalFlow_ECAdvection2D_t(this,s,nhat) result(exts) + !! Mirror boundary condition: sets extBoundary = interior state. + !! With the LLF Riemann flux, this gives sR = sL at the boundary, + !! so the Riemann flux reduces to the central flux (a.n)*s — no + !! dissipation. Use this BC when testing entropy conservation. + class(ECAdvection2D_t),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: nhat(1:2) + real(prec) :: exts(1:this%nvar) + + exts = s + if(.false.) exts(1) = exts(1)+nhat(1) ! suppress unused-dummy-argument warning + + endfunction hbc2d_NoNormalFlow_ECAdvection2D_t + +endmodule SELF_ECAdvection2D_t diff --git a/src/SELF_ECAdvection3D_t.f90 b/src/SELF_ECAdvection3D_t.f90 new file mode 100644 index 000000000..4247992c5 --- /dev/null +++ b/src/SELF_ECAdvection3D_t.f90 @@ -0,0 +1,123 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECAdvection3D_t + !! Entropy-conserving linear scalar advection equation in 3-D. + !! + !! Solves: + !! du/dt + u_a * du/dx + v_a * du/dy + w_a * du/dz = 0 + !! + !! Entropy function: eta(u) = u^2 / 2 + !! + !! EC two-point flux (arithmetic mean): + !! F^EC(uL, uR) = (u_a, v_a, w_a) * (uL+uR)/2 + !! + !! Surface Riemann flux (Local Lax-Friedrichs / Rusanov): + !! F_Riemann = 0.5 * (a.n) * (uL+uR) - 0.5 * |a| * (uR-uL) + !! where |a| = sqrt(u^2 + v^2 + w^2) is the maximum wave speed. + + use SELF_ECDGModel3D_t + use SELF_mesh + + implicit none + + type,extends(ECDGModel3D_t) :: ECAdvection3D_t + + real(prec) :: u ! x-component of advection velocity + real(prec) :: v ! y-component of advection velocity + real(prec) :: w ! z-component of advection velocity + + contains + + procedure :: entropy_func => entropy_func_ECAdvection3D_t + procedure :: twopointflux3d => twopointflux3d_ECAdvection3D_t + procedure :: riemannflux3d => riemannflux3d_ECAdvection3D_t + procedure :: hbc3d_NoNormalFlow => hbc3d_NoNormalFlow_ECAdvection3D_t + + endtype ECAdvection3D_t + +contains + + pure function entropy_func_ECAdvection3D_t(this,s) result(e) + !! Quadratic entropy: eta(u) = u^2 / 2 + class(ECAdvection3D_t),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec) :: e + + e = 0.5_prec*s(1)*s(1) + if(.false.) e = e+this%u ! suppress unused-dummy-argument warning + + endfunction entropy_func_ECAdvection3D_t + + pure function twopointflux3d_ECAdvection3D_t(this,sL,sR) result(flux) + !! Arithmetic-mean two-point flux for linear advection. + !! Entropy-conserving with respect to eta(u) = u^2/2. + class(ECAdvection3D_t),intent(in) :: this + real(prec),intent(in) :: sL(1:this%nvar) + real(prec),intent(in) :: sR(1:this%nvar) + real(prec) :: flux(1:this%nvar,1:3) + ! Local + real(prec) :: savg + + savg = 0.5_prec*(sL(1)+sR(1)) + flux(1,1) = this%u*savg + flux(1,2) = this%v*savg + flux(1,3) = this%w*savg + + endfunction twopointflux3d_ECAdvection3D_t + + pure function riemannflux3d_ECAdvection3D_t(this,sL,sR,dsdx,nhat) result(flux) + !! Local Lax-Friedrichs (Rusanov) Riemann flux for linear advection. + !! lambda_max = sqrt(u^2 + v^2 + w^2). + class(ECAdvection3D_t),intent(in) :: this + real(prec),intent(in) :: sL(1:this%nvar) + real(prec),intent(in) :: sR(1:this%nvar) + real(prec),intent(in) :: dsdx(1:this%nvar,1:3) + real(prec),intent(in) :: nhat(1:3) + real(prec) :: flux(1:this%nvar) + ! Local + real(prec) :: un,lam + + un = this%u*nhat(1)+this%v*nhat(2)+this%w*nhat(3) + lam = sqrt(this%u*this%u+this%v*this%v+this%w*this%w) + flux(1) = 0.5_prec*(un*(sL(1)+sR(1))-lam*(sR(1)-sL(1))) + if(.false.) flux(1) = flux(1)+dsdx(1,1) ! suppress unused-dummy-argument warning + + endfunction riemannflux3d_ECAdvection3D_t + + pure function hbc3d_NoNormalFlow_ECAdvection3D_t(this,s,nhat) result(exts) + !! Mirror boundary condition: sets extBoundary = interior state. + class(ECAdvection3D_t),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: nhat(1:3) + real(prec) :: exts(1:this%nvar) + + exts = s + if(.false.) exts(1) = exts(1)+nhat(1) ! suppress unused-dummy-argument warning + + endfunction hbc3d_NoNormalFlow_ECAdvection3D_t + +endmodule SELF_ECAdvection3D_t diff --git a/src/SELF_ECDGModel2D_t.f90 b/src/SELF_ECDGModel2D_t.f90 new file mode 100644 index 000000000..631e183df --- /dev/null +++ b/src/SELF_ECDGModel2D_t.f90 @@ -0,0 +1,234 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECDGModel2D_t + !! Entropy-conserving DGSEM base class for 2-D conservation laws. + !! + !! Extends DGModel2D_t by replacing the standard volume flux divergence with + !! the split-form (two-point) EC volume term following + !! Gassner, Winters, Kopriva (2016) + !! Winters, Kopriva, Gassner, Hindenlang (2021) + !! + !! Semidiscretization: + !! du/dt = source - (1/J) * [ 2 sum_n D_{n,i} F~^1(n,i,j) + !! + 2 sum_n D_{n,j} F~^2(n,i,j) + !! + M^{-1} B^T f_Riemann ] + !! + !! where F~^r is the Jacobian-weighted contravariant EC two-point flux in the + !! r-th computational direction (computed by TwoPointFluxMethod via twopointflux2d), + !! and f_Riemann is a standard Riemann solver flux at element faces (computed + !! by BoundaryFlux via riemannflux2d). + + use SELF_DGModel2D_t + use SELF_MappedTwoPointVector_2D + + implicit none + + type,extends(DGModel2D_t) :: ECDGModel2D_t + + type(MappedTwoPointVector2D) :: twoPointFlux + + contains + + procedure :: Init => Init_ECDGModel2D_t + procedure :: Free => Free_ECDGModel2D_t + + procedure :: twopointflux2d => twopointflux2d_ECDGModel2D_t + procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECDGModel2D_t + procedure :: CalculateTendency => CalculateTendency_ECDGModel2D_t + + endtype ECDGModel2D_t + +contains + + subroutine Init_ECDGModel2D_t(this,mesh,geometry) + implicit none + class(ECDGModel2D_t),intent(out) :: this + type(Mesh2D),intent(in),target :: mesh + type(SEMQuad),intent(in),target :: geometry + + ! Initialise all parent fields (solution, flux, source, ...) + call Init_DGModel2D_t(this,mesh,geometry) + + ! Additional two-point flux field + call this%twoPointFlux%Init(geometry%x%interp,this%nvar,mesh%nElem) + call this%twoPointFlux%AssociateGeometry(geometry) + + endsubroutine Init_ECDGModel2D_t + + subroutine Free_ECDGModel2D_t(this) + implicit none + class(ECDGModel2D_t),intent(inout) :: this + + call this%twoPointFlux%DissociateGeometry() + call this%twoPointFlux%Free() + call Free_DGModel2D_t(this) + + endsubroutine Free_ECDGModel2D_t + + pure function twopointflux2d_ECDGModel2D_t(this,sL,sR) result(flux) + !! Entropy-conserving two-point flux function. + !! + !! Returns the physical-space two-point EC flux for the state pair (sL, sR). + !! flux(ivar, d) is the d-th physical component (d=1: x, d=2: y) for variable ivar. + !! + !! Override this procedure in a concrete model to supply an EC flux. + !! The stub returns zero (no flux), which is safe but trivial. + implicit none + class(ECDGModel2D_t),intent(in) :: this + real(prec),intent(in) :: sL(1:this%nvar) + real(prec),intent(in) :: sR(1:this%nvar) + real(prec) :: flux(1:this%nvar,1:2) + + flux = 0.0_prec + if(.false.) flux(1,1) = flux(1,1)+sL(1)+sR(1) ! suppress unused-argument warning + + endfunction twopointflux2d_ECDGModel2D_t + + subroutine TwoPointFluxMethod_ECDGModel2D_t(this) + !! Computes pre-projected SCALAR contravariant two-point fluxes for all + !! node pairs and stores them in twoPointFlux%interior(n,i,j,iel,ivar,r). + !! + !! Following Trixi.jl for curved meshes, each computational direction r + !! uses the correct partner node AND the correct averaged metric Ja^r: + !! + !! r=1: pair (i,j)-(nn,j), Ja^1_avg = 0.5*(Ja^1(i,j) + Ja^1(nn,j)) + !! interior(nn,i,j,...,1) = sum_d Ja^1_avg(d) * F_EC_d(sL, sR) + !! + !! r=2: pair (i,j)-(i,nn), Ja^2_avg = 0.5*(Ja^2(i,j) + Ja^2(i,nn)) + !! interior(nn,i,j,...,2) = sum_d Ja^2_avg(d) * F_EC_d(sL, sR) + !! + !! The result is a SCALAR per (nn,i,j,iel,ivar,r) — no cross-contamination + !! between directions. + implicit none + class(ECDGModel2D_t),intent(inout) :: this + ! Local + integer :: nn,i,j,d,iEl,iVar + real(prec) :: sL(1:this%nvar),sR(1:this%nvar) + real(prec) :: Fphys(1:this%nvar,1:2) + real(prec) :: Fc + + do concurrent(nn=1:this%solution%N+1,i=1:this%solution%N+1, & + j=1:this%solution%N+1,iEl=1:this%mesh%nElem) + + sL = this%solution%interior(i,j,iEl,1:this%nvar) + + ! xi^1: pair (i,j)-(nn,j), project onto avg(Ja^1) + sR = this%solution%interior(nn,j,iEl,1:this%nvar) + Fphys = this%twopointflux2d(sL,sR) + do iVar = 1,this%nvar + Fc = 0.0_prec + do d = 1,2 + Fc = Fc+0.5_prec*( & + this%geometry%dsdx%interior(i,j,iEl,1,d,1)+ & + this%geometry%dsdx%interior(nn,j,iEl,1,d,1))* & + Fphys(iVar,d) + enddo + this%twoPointFlux%interior(nn,i,j,iEl,iVar,1) = Fc + enddo + + ! xi^2: pair (i,j)-(i,nn), project onto avg(Ja^2) + sR = this%solution%interior(i,nn,iEl,1:this%nvar) + Fphys = this%twopointflux2d(sL,sR) + do iVar = 1,this%nvar + Fc = 0.0_prec + do d = 1,2 + Fc = Fc+0.5_prec*( & + this%geometry%dsdx%interior(i,j,iEl,1,d,2)+ & + this%geometry%dsdx%interior(i,nn,iEl,1,d,2))* & + Fphys(iVar,d) + enddo + this%twoPointFlux%interior(nn,i,j,iEl,iVar,2) = Fc + enddo + + enddo + + endsubroutine TwoPointFluxMethod_ECDGModel2D_t + + subroutine CalculateTendency_ECDGModel2D_t(this) + !! Computes du/dt = source - EC-DG flux divergence. + !! + !! Volume term : MappedDivergence of twoPointFlux (EC split-form sum, /J) + !! Surface term: (1/J) * M^{-1} B^T f_Riemann (weak-form Riemann flux) + implicit none + class(ECDGModel2D_t),intent(inout) :: this + ! Local + integer :: i,j,iEl,iVar + real(prec) :: bM1,bM2,qwi,qwj,jac + + call this%solution%BoundaryInterp() + call this%solution%SideExchange(this%mesh) + + call this%PreTendency() + call this%SetBoundaryCondition() + + if(this%gradient_enabled) then + call this%CalculateSolutionGradient() + call this%SetGradientBoundaryCondition() + call this%solutionGradient%AverageSides() + endif + + call this%SourceMethod() + call this%BoundaryFlux() ! fills flux%boundaryNormal with Riemann fluxes + + ! EC volume divergence (includes J weighting via MappedDivergence) + call this%TwoPointFluxMethod() + call this%twoPointFlux%UpdateDevice() + call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior) + + ! Add weak-form DG surface term: (1/J) * M^{-1} B^T f_Riemann + ! Boundary side ordering: 1=South, 2=East, 3=North, 4=West + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar) + + bM1 = this%solution%interp%bMatrix(i,2) ! right/north boundary value at node i + bM2 = this%solution%interp%bMatrix(i,1) ! left/south boundary value at node i + qwi = this%solution%interp%qWeights(i) + qwj = this%solution%interp%qWeights(j) + jac = this%geometry%J%interior(i,j,iEl,1) + + ! East (side 2) and West (side 4) contributions — vary in xi^1 (i-direction) + this%fluxDivergence%interior(i,j,iEl,iVar) = & + this%fluxDivergence%interior(i,j,iEl,iVar)+ & + (bM1*this%flux%boundaryNormal(j,2,iEl,iVar)+ & + bM2*this%flux%boundaryNormal(j,4,iEl,iVar))/(qwi*jac)+ & + (this%solution%interp%bMatrix(j,2)*this%flux%boundaryNormal(i,3,iEl,iVar)+ & + this%solution%interp%bMatrix(j,1)*this%flux%boundaryNormal(i,1,iEl,iVar))/(qwj*jac) + + enddo + + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + iEl=1:this%mesh%nElem,iVar=1:this%solution%nVar) + + this%dSdt%interior(i,j,iEl,iVar) = & + this%source%interior(i,j,iEl,iVar)- & + this%fluxDivergence%interior(i,j,iEl,iVar) + + enddo + + endsubroutine CalculateTendency_ECDGModel2D_t + +endmodule SELF_ECDGModel2D_t diff --git a/src/SELF_ECDGModel3D_t.f90 b/src/SELF_ECDGModel3D_t.f90 new file mode 100644 index 000000000..5e004882a --- /dev/null +++ b/src/SELF_ECDGModel3D_t.f90 @@ -0,0 +1,221 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECDGModel3D_t + !! Entropy-conserving DGSEM base class for 3-D conservation laws. + !! See SELF_ECDGModel2D_t for design documentation. + + use SELF_DGModel3D_t + use SELF_MappedTwoPointVector_3D + + implicit none + + type,extends(DGModel3D_t) :: ECDGModel3D_t + + type(MappedTwoPointVector3D) :: twoPointFlux + + contains + + procedure :: Init => Init_ECDGModel3D_t + procedure :: Free => Free_ECDGModel3D_t + + procedure :: twopointflux3d => twopointflux3d_ECDGModel3D_t + procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECDGModel3D_t + procedure :: CalculateTendency => CalculateTendency_ECDGModel3D_t + + endtype ECDGModel3D_t + +contains + + subroutine Init_ECDGModel3D_t(this,mesh,geometry) + implicit none + class(ECDGModel3D_t),intent(out) :: this + type(Mesh3D),intent(in),target :: mesh + type(SEMHex),intent(in),target :: geometry + + call Init_DGModel3D_t(this,mesh,geometry) + + call this%twoPointFlux%Init(geometry%x%interp,this%nvar,mesh%nElem) + call this%twoPointFlux%AssociateGeometry(geometry) + + endsubroutine Init_ECDGModel3D_t + + subroutine Free_ECDGModel3D_t(this) + implicit none + class(ECDGModel3D_t),intent(inout) :: this + + call this%twoPointFlux%DissociateGeometry() + call this%twoPointFlux%Free() + call Free_DGModel3D_t(this) + + endsubroutine Free_ECDGModel3D_t + + pure function twopointflux3d_ECDGModel3D_t(this,sL,sR) result(flux) + !! Entropy-conserving two-point flux function (3-D). + !! flux(ivar, d) is the d-th physical component (d=1:x, d=2:y, d=3:z). + !! Override in concrete models. Stub returns zero. + implicit none + class(ECDGModel3D_t),intent(in) :: this + real(prec),intent(in) :: sL(1:this%nvar) + real(prec),intent(in) :: sR(1:this%nvar) + real(prec) :: flux(1:this%nvar,1:3) + + flux = 0.0_prec + if(.false.) flux(1,1) = flux(1,1)+sL(1)+sR(1) ! suppress unused-argument warning + + endfunction twopointflux3d_ECDGModel3D_t + + subroutine TwoPointFluxMethod_ECDGModel3D_t(this) + !! Computes pre-projected SCALAR contravariant two-point fluxes for all + !! node pairs, following Trixi.jl for curved meshes. Each direction r + !! uses the correct partner AND the correct averaged metric Ja^r. + implicit none + class(ECDGModel3D_t),intent(inout) :: this + ! Local + integer :: nn,i,j,k,d,iEl,iVar + real(prec) :: sL(1:this%nvar),sR(1:this%nvar) + real(prec) :: Fphys(1:this%nvar,1:3) + real(prec) :: Fc + + do concurrent(nn=1:this%solution%N+1,i=1:this%solution%N+1, & + j=1:this%solution%N+1,k=1:this%solution%N+1, & + iEl=1:this%mesh%nElem) + + sL = this%solution%interior(i,j,k,iEl,1:this%nvar) + + ! xi^1: pair (i,j,k)-(nn,j,k), project onto avg(Ja^1) + sR = this%solution%interior(nn,j,k,iEl,1:this%nvar) + Fphys = this%twopointflux3d(sL,sR) + do iVar = 1,this%nvar + Fc = 0.0_prec + do d = 1,3 + Fc = Fc+0.5_prec*( & + this%geometry%dsdx%interior(i,j,k,iEl,1,d,1)+ & + this%geometry%dsdx%interior(nn,j,k,iEl,1,d,1))* & + Fphys(iVar,d) + enddo + this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,1) = Fc + enddo + + ! xi^2: pair (i,j,k)-(i,nn,k), project onto avg(Ja^2) + sR = this%solution%interior(i,nn,k,iEl,1:this%nvar) + Fphys = this%twopointflux3d(sL,sR) + do iVar = 1,this%nvar + Fc = 0.0_prec + do d = 1,3 + Fc = Fc+0.5_prec*( & + this%geometry%dsdx%interior(i,j,k,iEl,1,d,2)+ & + this%geometry%dsdx%interior(i,nn,k,iEl,1,d,2))* & + Fphys(iVar,d) + enddo + this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,2) = Fc + enddo + + ! xi^3: pair (i,j,k)-(i,j,nn), project onto avg(Ja^3) + sR = this%solution%interior(i,j,nn,iEl,1:this%nvar) + Fphys = this%twopointflux3d(sL,sR) + do iVar = 1,this%nvar + Fc = 0.0_prec + do d = 1,3 + Fc = Fc+0.5_prec*( & + this%geometry%dsdx%interior(i,j,k,iEl,1,d,3)+ & + this%geometry%dsdx%interior(i,j,nn,iEl,1,d,3))* & + Fphys(iVar,d) + enddo + this%twoPointFlux%interior(nn,i,j,k,iEl,iVar,3) = Fc + enddo + + enddo + + endsubroutine TwoPointFluxMethod_ECDGModel3D_t + + subroutine CalculateTendency_ECDGModel3D_t(this) + implicit none + class(ECDGModel3D_t),intent(inout) :: this + ! Local + integer :: i,j,k,iEl,iVar + real(prec) :: bMi1,bMi2,bMj1,bMj2,bMk1,bMk2,qwi,qwj,qwk,jac + + call this%solution%BoundaryInterp() + call this%solution%SideExchange(this%mesh) + + call this%PreTendency() + call this%SetBoundaryCondition() + + if(this%gradient_enabled) then + call this%CalculateSolutionGradient() + call this%SetGradientBoundaryCondition() + call this%solutionGradient%AverageSides() + endif + + call this%SourceMethod() + call this%BoundaryFlux() + + call this%TwoPointFluxMethod() + call this%twoPointFlux%UpdateDevice() + call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior) + + ! Add (1/J) * M^{-1} B^T f_Riemann + ! 3D boundary side ordering: 1=Bottom, 2=South, 3=East, 4=North, 5=West, 6=Top + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + k=1:this%solution%N+1,iEl=1:this%mesh%nElem, & + iVar=1:this%solution%nVar) + + bMi1 = this%solution%interp%bMatrix(i,2) + bMi2 = this%solution%interp%bMatrix(i,1) + bMj1 = this%solution%interp%bMatrix(j,2) + bMj2 = this%solution%interp%bMatrix(j,1) + bMk1 = this%solution%interp%bMatrix(k,2) + bMk2 = this%solution%interp%bMatrix(k,1) + qwi = this%solution%interp%qWeights(i) + qwj = this%solution%interp%qWeights(j) + qwk = this%solution%interp%qWeights(k) + jac = this%geometry%J%interior(i,j,k,iEl,1) + + this%fluxDivergence%interior(i,j,k,iEl,iVar) = & + this%fluxDivergence%interior(i,j,k,iEl,iVar)+ & + (bMk1*this%flux%boundaryNormal(i,j,6,iEl,iVar)+ & ! top + bMk2*this%flux%boundaryNormal(i,j,1,iEl,iVar))/(qwk*jac)+ & ! bottom + (bMi1*this%flux%boundaryNormal(j,k,3,iEl,iVar)+ & ! east + bMi2*this%flux%boundaryNormal(j,k,5,iEl,iVar))/(qwi*jac)+ & ! west + (bMj1*this%flux%boundaryNormal(i,k,4,iEl,iVar)+ & ! north + bMj2*this%flux%boundaryNormal(i,k,2,iEl,iVar))/(qwj*jac) ! south + + enddo + + do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & + k=1:this%solution%N+1,iEl=1:this%mesh%nElem, & + iVar=1:this%solution%nVar) + + this%dSdt%interior(i,j,k,iEl,iVar) = & + this%source%interior(i,j,k,iEl,iVar)- & + this%fluxDivergence%interior(i,j,k,iEl,iVar) + + enddo + + endsubroutine CalculateTendency_ECDGModel3D_t + +endmodule SELF_ECDGModel3D_t diff --git a/src/SELF_Lagrange_t.f90 b/src/SELF_Lagrange_t.f90 index 6c0d85163..5469ff77c 100644 --- a/src/SELF_Lagrange_t.f90 +++ b/src/SELF_Lagrange_t.f90 @@ -97,6 +97,27 @@ module SELF_Lagrange_t !! The derivative matrix for mapping function nodal values to a nodal values of the derivative estimate. The dgMatrix is based !! on a weak form of the derivative. It must be used with bMatrix to account for boundary contributions in the weak form. + real(prec),pointer,contiguous,dimension(:,:) :: dSplitMatrix + !! The split-form derivative matrix D_split = D - 0.5*M^{-1}*B, where B is the SBP boundary operator. + !! + !! D_split is skew-symmetric under the M inner product: M*D_split + D_split^T*M = 0. Unlike D, it is + !! NOT an SBP operator (the weighted symmetric part is zero, not the boundary term B). This property + !! makes D_split ideal for the EC-DGSEM split-form volume integral: + !! + !! du/dt = -(2/J) sum_n D_split[n,i] * F_EC(u_i, u_n) + (1/J) M^{-1} B^T f_Riemann + !! + !! The skew-symmetry of the volume term guarantees it contributes zero to the entropy rate, so all + !! entropy change passes through the surface term. An entropy-dissipative Riemann solver then makes + !! the full scheme entropy-stable. + !! + !! Using D_split in the volume combined with plain f_Riemann on the surface is algebraically identical + !! to using D in the volume with the penalty (f_Riemann - f_local) on the surface (Trixi.jl convention). + !! + !! In SELF index convention (dSplitMatrix(ii,i) = D_split[i-1, ii-1]): + !! dSplitMatrix(ii,i) = dMatrix(ii,i) + !! - 0.5*(bMatrix(i,2)*bMatrix(ii,2) + !! - bMatrix(i,1)*bMatrix(ii,1)) / qWeights(i) + real(prec),pointer,contiguous,dimension(:,:) :: bMatrix !! The boundary interpolation matrix that is used to map a grid of nodal values at the control points to the element boundaries. @@ -149,6 +170,7 @@ subroutine Init_Lagrange_t(this,N,controlNodeType,M,targetNodeType) this%iMatrix(1:N+1,1:M+1), & this%dMatrix(1:N+1,1:N+1), & this%dgMatrix(1:N+1,1:N+1), & + this%dSplitMatrix(1:N+1,1:N+1), & this%bMatrix(1:N+1,1:2)) if(controlNodeType == GAUSS .or. controlNodeType == GAUSS_LOBATTO) then @@ -192,6 +214,21 @@ subroutine Init_Lagrange_t(this,N,controlNodeType,M,targetNodeType) this%bMatrix(1:N+1,1) = this%CalculateLagrangePolynomials(-1.0_prec) this%bMatrix(1:N+1,2) = this%CalculateLagrangePolynomials(1.0_prec) + ! D_split = D - 0.5*M^{-1}*B (skew-symmetric under the M inner product) + ! In SELF index convention: dSplitMatrix(ii,i) = dMatrix(ii,i) + ! - 0.5*(bMatrix(i,2)*bMatrix(ii,2) - bMatrix(i,1)*bMatrix(ii,1)) / qWeights(i) + block + integer :: i,ii + do i = 1,N+1 + do ii = 1,N+1 + this%dSplitMatrix(ii,i) = this%dMatrix(ii,i)- & + 0.5_prec*(this%bMatrix(i,2)*this%bMatrix(ii,2)- & + this%bMatrix(i,1)*this%bMatrix(ii,1))/ & + this%qWeights(i) + enddo + enddo + endblock + endsubroutine Init_Lagrange_t subroutine Free_Lagrange_t(this) @@ -207,6 +244,7 @@ subroutine Free_Lagrange_t(this) deallocate(this%iMatrix) deallocate(this%dMatrix) deallocate(this%dgMatrix) + deallocate(this%dSplitMatrix) deallocate(this%bMatrix) endsubroutine Free_Lagrange_t diff --git a/src/SELF_MappedTwoPointVector_2D_t.f90 b/src/SELF_MappedTwoPointVector_2D_t.f90 new file mode 100644 index 000000000..9002a631e --- /dev/null +++ b/src/SELF_MappedTwoPointVector_2D_t.f90 @@ -0,0 +1,130 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_MappedTwoPointVector_2D_t + !! Geometry-aware two-point vector and split-form divergence for 2-D. + !! + !! The MappedDivergence routine follows the construction described in + !! + !! Winters, Kopriva, Gassner, Hindenlang, + !! "Construction of Modern Robust Nodal Discontinuous Galerkin Spectral + !! Element Methods for the Compressible Navier-Stokes Equations", + !! Lecture Notes in Computational Science and Engineering, 2021. + !! + !! The key ingredient is that the contravariant two-point flux in the r-th + !! computational direction is formed by projecting the physical two-point + !! flux onto *averaged* contravariant basis vectors: + !! + !! F~^r_{(i,n),j} = sum_d (Ja^r_d(i,j) + Ja^r_d(n,j))/2 * f^d_{(i,n),j} + !! + !! where dsdx%interior(i,j,iEl,1,d,r) = J*a^r_d stores the scaled + !! contravariant basis vectors (metric terms times Jacobian). + !! The physical-space divergence at (i,j) is then + !! + !! (1/J_{i,j}) * 2 * sum_n [ D_{n,i} F~^1_{(i,n),j} + D_{n,j} F~^2_{i,(j,n)} ] + !! + !! The interior array stores the physical-space two-point flux + !! interior(n,i,j,iEl,iVar,d) = f^d_{(i or j, n)} where d is the physical + !! direction. Metric averaging is performed inside MappedDivergence. + + use SELF_Constants + use SELF_Lagrange + use SELF_Geometry_2D + use SELF_TwoPointVector_2D + use iso_c_binding + + implicit none + + type,extends(TwoPointVector2D),public :: MappedTwoPointVector2D_t + + logical :: geometry_associated = .false. + type(SEMQuad),pointer :: geometry => null() + + contains + + procedure,public :: AssociateGeometry => AssociateGeometry_MappedTwoPointVector2D_t + procedure,public :: DissociateGeometry => DissociateGeometry_MappedTwoPointVector2D_t + + generic,public :: MappedDivergence => MappedDivergence_MappedTwoPointVector2D_t + procedure,private :: MappedDivergence_MappedTwoPointVector2D_t + + endtype MappedTwoPointVector2D_t + +contains + + subroutine AssociateGeometry_MappedTwoPointVector2D_t(this,geometry) + implicit none + class(MappedTwoPointVector2D_t),intent(inout) :: this + type(SEMQuad),target,intent(in) :: geometry + + if(.not. associated(this%geometry)) then + this%geometry => geometry + this%geometry_associated = .true. + endif + + endsubroutine AssociateGeometry_MappedTwoPointVector2D_t + + subroutine DissociateGeometry_MappedTwoPointVector2D_t(this) + implicit none + class(MappedTwoPointVector2D_t),intent(inout) :: this + + if(associated(this%geometry)) then + this%geometry => null() + this%geometry_associated = .false. + endif + + endsubroutine DissociateGeometry_MappedTwoPointVector2D_t + + subroutine MappedDivergence_MappedTwoPointVector2D_t(this,df) + !! Computes the physical-space divergence of a 2-D split-form vector field + !! on a curvilinear mesh. + !! + !! Convention (following Trixi.jl for curved meshes): + !! interior(n,i,j,iEl,iVar,r) holds the pre-projected SCALAR contravariant + !! two-point flux for the r-th computational direction: + !! + !! interior(n,i,j,iEl,iVar,1) = avg(Ja^1) . F_EC(s(i,j), s(n,j)) + !! interior(n,i,j,iEl,iVar,2) = avg(Ja^2) . F_EC(s(i,j), s(i,n)) + !! + !! The metric averaging and direction-correct pairing are the caller's + !! responsibility (e.g. TwoPointFluxMethod in ECDGModel). + !! MappedDivergence applies the reference-element split-form sum and + !! divides by J. + implicit none + class(MappedTwoPointVector2D_t),intent(in) :: this + real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%nElem,1:this%nVar) + ! Local + integer :: i,j,iEl,iVar + + call this%Divergence(df) + + do concurrent(i=1:this%N+1,j=1:this%N+1,iEl=1:this%nElem,iVar=1:this%nVar) + df(i,j,iEl,iVar) = df(i,j,iEl,iVar)/this%geometry%J%interior(i,j,iEl,1) + enddo + + endsubroutine MappedDivergence_MappedTwoPointVector2D_t + +endmodule SELF_MappedTwoPointVector_2D_t diff --git a/src/SELF_MappedTwoPointVector_3D_t.f90 b/src/SELF_MappedTwoPointVector_3D_t.f90 new file mode 100644 index 000000000..10c7b56e6 --- /dev/null +++ b/src/SELF_MappedTwoPointVector_3D_t.f90 @@ -0,0 +1,122 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_MappedTwoPointVector_3D_t + !! Geometry-aware two-point vector and split-form divergence for 3-D. + !! + !! The MappedDivergence routine follows the construction described in + !! + !! Winters, Kopriva, Gassner, Hindenlang, + !! "Construction of Modern Robust Nodal Discontinuous Galerkin Spectral + !! Element Methods for the Compressible Navier-Stokes Equations", + !! Lecture Notes in Computational Science and Engineering, 2021. + !! + !! The contravariant two-point flux in the r-th computational direction is + !! + !! F~^r_{(i,n),j,k} = sum_d (Ja^r_d(i,j,k) + Ja^r_d(n,j,k))/2 * f^d_{(i,n),j,k} + !! + !! and similarly for the xi^2 and xi^3 sums (averaging (i,j,k)-(i,n,k) and + !! (i,j,k)-(i,j,n) respectively). The physical divergence is + !! + !! (1/J_{i,j,k}) * 2 * sum_n [ D_{n,i} F~^1 + D_{n,j} F~^2 + D_{n,k} F~^3 ] + !! + !! interior(n,i,j,k,iEl,iVar,d) stores the physical-space two-point flux + !! f^d in the d-th physical direction. + + use SELF_Constants + use SELF_Lagrange + use SELF_Geometry_3D + use SELF_TwoPointVector_3D + use iso_c_binding + + implicit none + + type,extends(TwoPointVector3D),public :: MappedTwoPointVector3D_t + + logical :: geometry_associated = .false. + type(SEMHex),pointer :: geometry => null() + + contains + + procedure,public :: AssociateGeometry => AssociateGeometry_MappedTwoPointVector3D_t + procedure,public :: DissociateGeometry => DissociateGeometry_MappedTwoPointVector3D_t + + generic,public :: MappedDivergence => MappedDivergence_MappedTwoPointVector3D_t + procedure,private :: MappedDivergence_MappedTwoPointVector3D_t + + endtype MappedTwoPointVector3D_t + +contains + + subroutine AssociateGeometry_MappedTwoPointVector3D_t(this,geometry) + implicit none + class(MappedTwoPointVector3D_t),intent(inout) :: this + type(SEMHex),target,intent(in) :: geometry + + if(.not. associated(this%geometry)) then + this%geometry => geometry + this%geometry_associated = .true. + endif + + endsubroutine AssociateGeometry_MappedTwoPointVector3D_t + + subroutine DissociateGeometry_MappedTwoPointVector3D_t(this) + implicit none + class(MappedTwoPointVector3D_t),intent(inout) :: this + + if(associated(this%geometry)) then + this%geometry => null() + this%geometry_associated = .false. + endif + + endsubroutine DissociateGeometry_MappedTwoPointVector3D_t + + subroutine MappedDivergence_MappedTwoPointVector3D_t(this,df) + !! Computes the physical-space divergence of a 3-D split-form vector field + !! on a curvilinear mesh. + !! + !! Convention (following Trixi.jl for curved meshes): + !! interior(n,i,j,k,iEl,iVar,r) holds the pre-projected SCALAR contravariant + !! two-point flux for the r-th computational direction. + !! MappedDivergence applies the reference-element split-form sum and + !! divides by J. + implicit none + class(MappedTwoPointVector3D_t),intent(in) :: this + real(prec),intent(out) :: & + df(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nElem,1:this%nVar) + ! Local + integer :: i,j,k,iEl,iVar + + call this%Divergence(df) + + do concurrent(i=1:this%N+1,j=1:this%N+1,k=1:this%N+1, & + iEl=1:this%nElem,iVar=1:this%nVar) + df(i,j,k,iEl,iVar) = df(i,j,k,iEl,iVar)/this%geometry%J%interior(i,j,k,iEl,1) + enddo + + endsubroutine MappedDivergence_MappedTwoPointVector3D_t + +endmodule SELF_MappedTwoPointVector_3D_t diff --git a/src/SELF_TwoPointVector_2D_t.f90 b/src/SELF_TwoPointVector_2D_t.f90 new file mode 100644 index 000000000..7276f723b --- /dev/null +++ b/src/SELF_TwoPointVector_2D_t.f90 @@ -0,0 +1,173 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_TwoPointVector_2D_t + + use SELF_Constants + use SELF_Lagrange + use SELF_Metadata + use FEQParse + use SELF_Data + + use iso_c_binding + + implicit none + + type,extends(SELF_DataObj),public :: TwoPointVector2D_t + !! A two-point vector field for use in split-form DGSEM in 2-D. + !! + !! Memory layout: interior(n, i, j, nEl, nVar, 1:2) + !! dim 1 (n) : two-point index; looping over n is equivalent to looping + !! over n in the split-form divergence sum (0:N in 0-based) + !! dim 2 (i) : first computational coordinate direction (xi^1) + !! dim 3 (j) : second computational coordinate direction (xi^2) + !! dim 4 : element index + !! dim 5 : variable index + !! dim 6 (idir): vector direction (1 or 2) + + real(prec),pointer,contiguous,dimension(:,:,:,:,:,:) :: interior + + contains + + procedure,public :: Init => Init_TwoPointVector2D_t + procedure,public :: Free => Free_TwoPointVector2D_t + + procedure,public :: UpdateHost => UpdateHost_TwoPointVector2D_t + procedure,public :: UpdateDevice => UpdateDevice_TwoPointVector2D_t + + generic,public :: Divergence => Divergence_TwoPointVector2D_t + procedure,private :: Divergence_TwoPointVector2D_t + + endtype TwoPointVector2D_t + +contains + + subroutine Init_TwoPointVector2D_t(this,interp,nVar,nElem) + !! Allocate the interior array for a 2-D two-point vector field. + !! The interior array has rank 6 with layout (n,i,j,nEl,nVar,idir). + !! + !! Requires Gauss-Lobatto quadrature nodes (controlNodeType=GAUSS_LOBATTO). + !! With GLL nodes D_split = (D-D^T)/2 is exactly skew-symmetric, which is + !! the standard operator used in the entropy-conserving split-form DGSEM + !! (Gassner, Winters, Kopriva 2016). + implicit none + class(TwoPointVector2D_t),intent(out) :: this + type(Lagrange),target,intent(in) :: interp + integer,intent(in) :: nVar + integer,intent(in) :: nElem + ! Local + integer :: i + + if(interp%controlNodeType /= GAUSS_LOBATTO) then + print*,__FILE__//" : TwoPointVector2D requires Gauss-Lobatto quadrature nodes." + stop 1 + endif + + this%interp => interp + this%nVar = nVar + this%nElem = nElem + this%N = interp%N + this%M = interp%M + + allocate(this%interior(1:interp%N+1,1:interp%N+1,1:interp%N+1, & + 1:nElem,1:nVar,1:2)) + + allocate(this%meta(1:nVar)) + allocate(this%eqn(1:2*nVar)) + + ! Initialize equation parser to prevent segmentation faults with amdflang + ! when the parser functions are not allocated (see SELF_Vector_2D_t.f90) + do i = 1,2*nVar + this%eqn(i) = EquationParser('f=0',(/'x','y','z','t'/)) + enddo + + this%interior = 0.0_prec + + endsubroutine Init_TwoPointVector2D_t + + subroutine Free_TwoPointVector2D_t(this) + implicit none + class(TwoPointVector2D_t),intent(inout) :: this + + this%interp => null() + this%nVar = 0 + this%nElem = 0 + + deallocate(this%interior) + deallocate(this%meta) + deallocate(this%eqn) + + endsubroutine Free_TwoPointVector2D_t + + subroutine UpdateHost_TwoPointVector2D_t(this) + implicit none + class(TwoPointVector2D_t),intent(inout) :: this + if(.false.) this%N = this%N ! CPU stub; suppress unused-dummy-argument warning + endsubroutine UpdateHost_TwoPointVector2D_t + + subroutine UpdateDevice_TwoPointVector2D_t(this) + implicit none + class(TwoPointVector2D_t),intent(inout) :: this + if(.false.) this%N = this%N ! CPU stub; suppress unused-dummy-argument warning + endsubroutine UpdateDevice_TwoPointVector2D_t + + subroutine Divergence_TwoPointVector2D_t(this,df) + !! Computes the split-form (two-point) divergence of a 2-D vector field + !! in the reference element (computational coordinates). + !! + !! The split-form divergence at node (i,j) is + !! + !! (nabla.F)_{i,j} = 2 sum_n [ D_{n,i} F^1(n,i,j) + D_{n,j} F^2(n,i,j) ] + !! + !! where D is the standard derivative matrix (dMatrix) and + !! F^idir(n,i,j,...) = interior(n,i,j,iEl,iVar,idir) stores the two-point + !! flux between nodes i (or j) and n in the idir-th coordinate direction. + !! + !! The interior array is assumed to hold contravariant (J-scaled) two-point + !! fluxes. To obtain the physical divergence, divide the result by the + !! element Jacobian J(i,j,iEl). + implicit none + class(TwoPointVector2D_t),intent(in) :: this + real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%nElem,1:this%nVar) + ! Local + integer :: i,j,nn,iEl,iVar + real(prec) :: dfLoc + + do concurrent(i=1:this%N+1,j=1:this%N+1,iEl=1:this%nElem,iVar=1:this%nVar) + + dfLoc = 0.0_prec + do nn = 1,this%N+1 + dfLoc = dfLoc+ & + this%interp%dSplitMatrix(nn,i)*this%interior(nn,i,j,iEl,iVar,1)+ & + this%interp%dSplitMatrix(nn,j)*this%interior(nn,i,j,iEl,iVar,2) + enddo + df(i,j,iEl,iVar) = 2.0_prec*dfLoc + + enddo + + endsubroutine Divergence_TwoPointVector2D_t + +endmodule SELF_TwoPointVector_2D_t diff --git a/src/SELF_TwoPointVector_3D_t.f90 b/src/SELF_TwoPointVector_3D_t.f90 new file mode 100644 index 000000000..35845accd --- /dev/null +++ b/src/SELF_TwoPointVector_3D_t.f90 @@ -0,0 +1,175 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_TwoPointVector_3D_t + + use SELF_Constants + use SELF_Lagrange + use SELF_Metadata + use FEQParse + use SELF_Data + + use iso_c_binding + + implicit none + + type,extends(SELF_DataObj),public :: TwoPointVector3D_t + !! A two-point vector field for use in split-form DGSEM in 3-D. + !! + !! Memory layout: interior(n, i, j, k, nEl, nVar, 1:3) + !! dim 1 (n) : two-point index; looping over n is equivalent to looping + !! over n in the split-form divergence sum (0:N in 0-based) + !! dim 2 (i) : first computational coordinate direction (xi^1) + !! dim 3 (j) : second computational coordinate direction (xi^2) + !! dim 4 (k) : third computational coordinate direction (xi^3) + !! dim 5 : element index + !! dim 6 : variable index + !! dim 7 (idir): vector direction (1, 2, or 3) + + real(prec),pointer,contiguous,dimension(:,:,:,:,:,:,:) :: interior + + contains + + procedure,public :: Init => Init_TwoPointVector3D_t + procedure,public :: Free => Free_TwoPointVector3D_t + + procedure,public :: UpdateHost => UpdateHost_TwoPointVector3D_t + procedure,public :: UpdateDevice => UpdateDevice_TwoPointVector3D_t + + generic,public :: Divergence => Divergence_TwoPointVector3D_t + procedure,private :: Divergence_TwoPointVector3D_t + + endtype TwoPointVector3D_t + +contains + + subroutine Init_TwoPointVector3D_t(this,interp,nVar,nElem) + !! Allocate the interior array for a 3-D two-point vector field. + !! The interior array has rank 7 with layout (n,i,j,k,nEl,nVar,idir). + !! + !! Requires Gauss-Lobatto quadrature nodes (controlNodeType=GAUSS_LOBATTO). + implicit none + class(TwoPointVector3D_t),intent(out) :: this + type(Lagrange),target,intent(in) :: interp + integer,intent(in) :: nVar + integer,intent(in) :: nElem + ! Local + integer :: i + + if(interp%controlNodeType /= GAUSS_LOBATTO) then + print*,__FILE__//" : TwoPointVector3D requires Gauss-Lobatto quadrature nodes." + stop 1 + endif + + this%interp => interp + this%nVar = nVar + this%nElem = nElem + this%N = interp%N + this%M = interp%M + + allocate(this%interior(1:interp%N+1,1:interp%N+1,1:interp%N+1,1:interp%N+1, & + 1:nElem,1:nVar,1:3)) + + allocate(this%meta(1:nVar)) + allocate(this%eqn(1:3*nVar)) + + ! Initialize equation parser to prevent segmentation faults with amdflang + ! when the parser functions are not allocated (see SELF_Vector_2D_t.f90) + do i = 1,3*nVar + this%eqn(i) = EquationParser('f=0',(/'x','y','z','t'/)) + enddo + + this%interior = 0.0_prec + + endsubroutine Init_TwoPointVector3D_t + + subroutine Free_TwoPointVector3D_t(this) + implicit none + class(TwoPointVector3D_t),intent(inout) :: this + + this%interp => null() + this%nVar = 0 + this%nElem = 0 + + deallocate(this%interior) + deallocate(this%meta) + deallocate(this%eqn) + + endsubroutine Free_TwoPointVector3D_t + + subroutine UpdateHost_TwoPointVector3D_t(this) + implicit none + class(TwoPointVector3D_t),intent(inout) :: this + if(.false.) this%N = this%N ! CPU stub; suppress unused-dummy-argument warning + endsubroutine UpdateHost_TwoPointVector3D_t + + subroutine UpdateDevice_TwoPointVector3D_t(this) + implicit none + class(TwoPointVector3D_t),intent(inout) :: this + if(.false.) this%N = this%N ! CPU stub; suppress unused-dummy-argument warning + endsubroutine UpdateDevice_TwoPointVector3D_t + + subroutine Divergence_TwoPointVector3D_t(this,df) + !! Computes the split-form (two-point) divergence of a 3-D vector field + !! in the reference element (computational coordinates). + !! + !! The split-form divergence at node (i,j,k) is + !! + !! (nabla.F)_{i,j,k} = 2 sum_n [ D_{n,i} F^1(n,i,j,k) + !! + D_{n,j} F^2(n,i,j,k) + !! + D_{n,k} F^3(n,i,j,k) ] + !! + !! where D is the standard derivative matrix (dMatrix) and + !! F^idir(n,i,j,k,...) = interior(n,i,j,k,iEl,iVar,idir) stores the + !! two-point flux between nodes i (j or k) and n in the idir-th direction. + !! + !! The interior array is assumed to hold contravariant (J-scaled) two-point + !! fluxes. To obtain the physical divergence, divide the result by the + !! element Jacobian J(i,j,k,iEl). + implicit none + class(TwoPointVector3D_t),intent(in) :: this + real(prec),intent(out) :: df(1:this%N+1,1:this%N+1,1:this%N+1,1:this%nElem,1:this%nVar) + ! Local + integer :: i,j,k,nn,iEl,iVar + real(prec) :: dfLoc + + do concurrent(i=1:this%N+1,j=1:this%N+1,k=1:this%N+1, & + iEl=1:this%nElem,iVar=1:this%nVar) + + dfLoc = 0.0_prec + do nn = 1,this%N+1 + dfLoc = dfLoc+ & + this%interp%dSplitMatrix(nn,i)*this%interior(nn,i,j,k,iEl,iVar,1)+ & + this%interp%dSplitMatrix(nn,j)*this%interior(nn,i,j,k,iEl,iVar,2)+ & + this%interp%dSplitMatrix(nn,k)*this%interior(nn,i,j,k,iEl,iVar,3) + enddo + df(i,j,k,iEl,iVar) = 2.0_prec*dfLoc + + enddo + + endsubroutine Divergence_TwoPointVector3D_t + +endmodule SELF_TwoPointVector_3D_t diff --git a/src/cpu/SELF_ECAdvection2D.f90 b/src/cpu/SELF_ECAdvection2D.f90 new file mode 100644 index 000000000..6551d6afb --- /dev/null +++ b/src/cpu/SELF_ECAdvection2D.f90 @@ -0,0 +1,36 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECAdvection2D + + use SELF_ECAdvection2D_t + + implicit none + + type,extends(ECAdvection2D_t),public :: ECAdvection2D + endtype ECAdvection2D + +endmodule SELF_ECAdvection2D diff --git a/src/cpu/SELF_ECAdvection3D.f90 b/src/cpu/SELF_ECAdvection3D.f90 new file mode 100644 index 000000000..919f5561a --- /dev/null +++ b/src/cpu/SELF_ECAdvection3D.f90 @@ -0,0 +1,36 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECAdvection3D + + use SELF_ECAdvection3D_t + + implicit none + + type,extends(ECAdvection3D_t),public :: ECAdvection3D + endtype ECAdvection3D + +endmodule SELF_ECAdvection3D diff --git a/src/cpu/SELF_ECDGModel2D.f90 b/src/cpu/SELF_ECDGModel2D.f90 new file mode 100644 index 000000000..4d173ea0e --- /dev/null +++ b/src/cpu/SELF_ECDGModel2D.f90 @@ -0,0 +1,36 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECDGModel2D + + use SELF_ECDGModel2D_t + + implicit none + + type,extends(ECDGModel2D_t),public :: ECDGModel2D + endtype ECDGModel2D + +endmodule SELF_ECDGModel2D diff --git a/src/cpu/SELF_ECDGModel3D.f90 b/src/cpu/SELF_ECDGModel3D.f90 new file mode 100644 index 000000000..a2c84ee51 --- /dev/null +++ b/src/cpu/SELF_ECDGModel3D.f90 @@ -0,0 +1,36 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECDGModel3D + + use SELF_ECDGModel3D_t + + implicit none + + type,extends(ECDGModel3D_t),public :: ECDGModel3D + endtype ECDGModel3D + +endmodule SELF_ECDGModel3D diff --git a/src/cpu/SELF_MappedTwoPointVector_2D.f90 b/src/cpu/SELF_MappedTwoPointVector_2D.f90 new file mode 100644 index 000000000..a95eba870 --- /dev/null +++ b/src/cpu/SELF_MappedTwoPointVector_2D.f90 @@ -0,0 +1,36 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_MappedTwoPointVector_2D + + use SELF_MappedTwoPointVector_2D_t + + implicit none + + type,extends(MappedTwoPointVector2D_t),public :: MappedTwoPointVector2D + endtype MappedTwoPointVector2D + +endmodule SELF_MappedTwoPointVector_2D diff --git a/src/cpu/SELF_MappedTwoPointVector_3D.f90 b/src/cpu/SELF_MappedTwoPointVector_3D.f90 new file mode 100644 index 000000000..c6b9808c8 --- /dev/null +++ b/src/cpu/SELF_MappedTwoPointVector_3D.f90 @@ -0,0 +1,36 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_MappedTwoPointVector_3D + + use SELF_MappedTwoPointVector_3D_t + + implicit none + + type,extends(MappedTwoPointVector3D_t),public :: MappedTwoPointVector3D + endtype MappedTwoPointVector3D + +endmodule SELF_MappedTwoPointVector_3D diff --git a/src/cpu/SELF_TwoPointVector_2D.f90 b/src/cpu/SELF_TwoPointVector_2D.f90 new file mode 100644 index 000000000..511123a40 --- /dev/null +++ b/src/cpu/SELF_TwoPointVector_2D.f90 @@ -0,0 +1,37 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_TwoPointVector_2D + + use SELF_TwoPointVector_2D_t + + implicit none + + type,extends(TwoPointVector2D_t),public :: TwoPointVector2D + character(3) :: backend = "cpu" + endtype TwoPointVector2D + +endmodule SELF_TwoPointVector_2D diff --git a/src/cpu/SELF_TwoPointVector_3D.f90 b/src/cpu/SELF_TwoPointVector_3D.f90 new file mode 100644 index 000000000..e99b9bed2 --- /dev/null +++ b/src/cpu/SELF_TwoPointVector_3D.f90 @@ -0,0 +1,37 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_TwoPointVector_3D + + use SELF_TwoPointVector_3D_t + + implicit none + + type,extends(TwoPointVector3D_t),public :: TwoPointVector3D + character(3) :: backend = "cpu" + endtype TwoPointVector3D + +endmodule SELF_TwoPointVector_3D diff --git a/src/gpu/SELF_DGModel2D.f90 b/src/gpu/SELF_DGModel2D.f90 index 0d5c00042..aa4534912 100644 --- a/src/gpu/SELF_DGModel2D.f90 +++ b/src/gpu/SELF_DGModel2D.f90 @@ -166,7 +166,9 @@ subroutine CalculateEntropy_DGModel2D(this) do i = 1,this%solution%interp%N+1 jac = abs(this%geometry%J%interior(i,j,iel,1)) s = this%solution%interior(i,j,iel,1:this%nvar) - e = e+this%entropy_func(s)*jac + e = e+this%entropy_func(s)*jac* & + this%solution%interp%qWeights(i)* & + this%solution%interp%qWeights(j) enddo enddo enddo diff --git a/src/gpu/SELF_DGModel3D.f90 b/src/gpu/SELF_DGModel3D.f90 index 2fe9378f5..908f16b09 100644 --- a/src/gpu/SELF_DGModel3D.f90 +++ b/src/gpu/SELF_DGModel3D.f90 @@ -171,7 +171,10 @@ subroutine CalculateEntropy_DGModel3D(this) do i = 1,this%solution%interp%N+1 jac = abs(this%geometry%J%interior(i,j,k,iel,1)) s = this%solution%interior(i,j,k,iel,1:this%nvar) - e = e+this%entropy_func(s)*jac + e = e+this%entropy_func(s)*jac* & + this%solution%interp%qWeights(i)* & + this%solution%interp%qWeights(j)* & + this%solution%interp%qWeights(k) enddo enddo enddo diff --git a/src/gpu/SELF_ECAdvection2D.cpp b/src/gpu/SELF_ECAdvection2D.cpp new file mode 100644 index 000000000..f821cd8e7 --- /dev/null +++ b/src/gpu/SELF_ECAdvection2D.cpp @@ -0,0 +1,150 @@ +#include "SELF_GPU_Macros.h" +#include + +// ============================================================ +// Mirror boundary condition: extBoundary = boundary at domain faces +// ============================================================ + +__global__ void setboundarycondition_ecadvection2d_gpukernel( + real *extBoundary, real *boundary, int *sideInfo, int N, int nel, int nvar) +{ + uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; + uint32_t ndof = (N+1)*4*nel; + + if (idof < ndof) { + uint32_t i = idof % (N+1); + uint32_t s1 = (idof/(N+1)) % 4; + uint32_t e1 = idof/(N+1)/4; + uint32_t e2 = sideInfo[INDEX3(2,s1,e1,5,4)]; + if (e2 == 0) { + uint32_t ivar = blockIdx.y; + extBoundary[SCB_2D_INDEX(i,s1,e1,ivar,N,nel)] = + boundary[SCB_2D_INDEX(i,s1,e1,ivar,N,nel)]; + } + } +} + +extern "C" +{ + void setboundarycondition_ecadvection2d_gpu( + real *extBoundary, real *boundary, int *sideInfo, int N, int nel, int nvar) + { + int ndof = (N+1)*4*nel; + int threads_per_block = 256; + int nblocks_x = ndof/threads_per_block + 1; + setboundarycondition_ecadvection2d_gpukernel<<>>(extBoundary, boundary, sideInfo, N, nel, nvar); + } +} + +// ============================================================ +// LLF boundary flux for linear advection +// flux = 0.5*(un*(sL+sR) - lambda*(sR-sL)) * nScale +// un = u*nx + v*ny, lambda = sqrt(u^2+v^2) +// ============================================================ + +__global__ void boundaryflux_ecadvection2d_gpukernel( + real *fb, real *fextb, real *nhat, real *nscale, real *flux, + real u, real v, real lam, int N, int nel, int nvar) +{ + uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; + uint32_t ndof = (N+1)*4*nel; + + if (idof < ndof) { + uint32_t i = idof % (N+1); + uint32_t j = (idof/(N+1)) % 4; + uint32_t iel = idof/(N+1)/4; + uint32_t ivar = blockIdx.y; + + real nx = nhat[VEB_2D_INDEX(i,j,iel,0,0,N,nel,1)]; + real ny = nhat[VEB_2D_INDEX(i,j,iel,0,1,N,nel,1)]; + real un = u*nx + v*ny; + real nmag = nscale[SCB_2D_INDEX(i,j,iel,0,N,nel)]; + + real sL = fb[idof + ivar*ndof]; + real sR = fextb[idof + ivar*ndof]; + + flux[idof + ivar*ndof] = 0.5*(un*(sL+sR) - lam*(sR-sL)) * nmag; + } +} + +extern "C" +{ + void boundaryflux_ecadvection2d_gpu( + real *fb, real *fextb, real *nhat, real *nscale, real *flux, + real u, real v, int N, int nel, int nvar) + { + real lam = sqrt(u*u + v*v); + int ndof = (N+1)*4*nel; + int threads_per_block = 256; + int nblocks_x = ndof/threads_per_block + 1; + boundaryflux_ecadvection2d_gpukernel<<>>(fb, fextb, nhat, nscale, flux, u, v, lam, N, nel, nvar); + } +} + +// ============================================================ +// TwoPointFluxMethod for EC linear advection on a curvilinear mesh +// +// Computes contravariant two-point fluxes: +// Fc^r(nn,i,j) = avg(Ja^r . a) * (s(node_L) + s(node_R)) / 2 +// +// where a = (u,v) is the constant advection velocity and the metric +// averaging uses the correct partner per direction. +// +// f : output, TPV_2D_INDEX layout +// s : solution interior, SC_2D_INDEX layout +// dsdx: metric terms, TE_2D_INDEX layout (nVar=1 for geometry) +// ============================================================ + +template +__global__ void __launch_bounds__(256) twopointfluxmethod_ecadvection2d_gpukernel( + real *f, real *s, real *dsdx, real u, real v, int nq, int N, int nvar) +{ + uint32_t iq = threadIdx.x; + if (iq < nq) { + uint32_t iel = blockIdx.x; + uint32_t nel = gridDim.x; + uint32_t ivar = blockIdx.y; + uint32_t i = iq % (N+1); + uint32_t j = iq / (N+1); + + real s_ij = s[iq + nq*(iel + nel*ivar)]; + int nq3 = nq*(N+1); // stride between idir slices in TPV layout + + for (int nn = 0; nn < N+1; nn++) { + uint32_t iq_nn_j = nn + (N+1)*j; // node (nn,j) for xi^1 + uint32_t iq_i_nn = i + (N+1)*nn; // node (i,nn) for xi^2 + + // xi^1: pair (i,j)-(nn,j), project onto avg(Ja^1) + real s_nn_j = s[iq_nn_j + nq*(iel + nel*ivar)]; + real savg1 = 0.5*(s_ij + s_nn_j); + real un_contra1 = 0.5*(dsdx[iq + nq*(iel + nel*0)] + dsdx[iq_nn_j + nq*(iel + nel*0)])*u + + 0.5*(dsdx[iq + nq*(iel + nel*1)] + dsdx[iq_nn_j + nq*(iel + nel*1)])*v; + f[nn + (N+1)*iq + nq3*(iel + nel*ivar)] = un_contra1 * savg1; + + // xi^2: pair (i,j)-(i,nn), project onto avg(Ja^2) + real s_i_nn = s[iq_i_nn + nq*(iel + nel*ivar)]; + real savg2 = 0.5*(s_ij + s_i_nn); + real un_contra2 = 0.5*(dsdx[iq + nq*(iel + nel*2)] + dsdx[iq_i_nn + nq*(iel + nel*2)])*u + + 0.5*(dsdx[iq + nq*(iel + nel*3)] + dsdx[iq_i_nn + nq*(iel + nel*3)])*v; + f[nn + (N+1)*iq + nq3*(iel + nel*(ivar + nvar))] = un_contra2 * savg2; + } + } +} + +extern "C" +{ + void twopointfluxmethod_ecadvection2d_gpu( + real *f, real *s, real *dsdx, real u, real v, int N, int nvar, int nel) + { + int nq = (N+1)*(N+1); + if (N <= 7) { + twopointfluxmethod_ecadvection2d_gpukernel<64><<>>(f, s, dsdx, u, v, nq, N, nvar); + } else { + twopointfluxmethod_ecadvection2d_gpukernel<256><<>>(f, s, dsdx, u, v, nq, N, nvar); + } + } +} diff --git a/src/gpu/SELF_ECAdvection2D.f90 b/src/gpu/SELF_ECAdvection2D.f90 new file mode 100644 index 000000000..515ca9feb --- /dev/null +++ b/src/gpu/SELF_ECAdvection2D.f90 @@ -0,0 +1,141 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECAdvection2D + + use SELF_ECAdvection2D_t + use SELF_GPU + use SELF_GPUInterfaces + use iso_c_binding + + implicit none + + type,extends(ECAdvection2D_t),public :: ECAdvection2D + + contains + + procedure :: SetBoundaryCondition => SetBoundaryCondition_ECAdvection2D + procedure :: BoundaryFlux => BoundaryFlux_ECAdvection2D + procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECAdvection2D + procedure :: SourceMethod => SourceMethod_ECAdvection2D + + endtype ECAdvection2D + + interface + subroutine setboundarycondition_ecadvection2d_gpu(extboundary,boundary,sideinfo,N,nel,nvar) & + bind(c,name="setboundarycondition_ecadvection2d_gpu") + use iso_c_binding + type(c_ptr),value :: extboundary,boundary,sideinfo + integer(c_int),value :: N,nel,nvar + endsubroutine setboundarycondition_ecadvection2d_gpu + endinterface + + interface + subroutine boundaryflux_ecadvection2d_gpu(fb,fextb,nhat,nscale,flux,u,v,N,nel,nvar) & + bind(c,name="boundaryflux_ecadvection2d_gpu") + use iso_c_binding + use SELF_Constants + type(c_ptr),value :: fb,fextb,nhat,nscale,flux + real(c_prec),value :: u,v + integer(c_int),value :: N,nel,nvar + endsubroutine boundaryflux_ecadvection2d_gpu + endinterface + + interface + subroutine twopointfluxmethod_ecadvection2d_gpu(f,s,dsdx,u,v,N,nvar,nel) & + bind(c,name="twopointfluxmethod_ecadvection2d_gpu") + use iso_c_binding + use SELF_Constants + type(c_ptr),value :: f,s,dsdx + real(c_prec),value :: u,v + integer(c_int),value :: N,nvar,nel + endsubroutine twopointfluxmethod_ecadvection2d_gpu + endinterface + +contains + + subroutine SetBoundaryCondition_ECAdvection2D(this) + !! Mirror BC on GPU: extBoundary = boundary at all domain faces. + implicit none + class(ECAdvection2D),intent(inout) :: this + + call setboundarycondition_ecadvection2d_gpu( & + this%solution%extboundary_gpu, & + this%solution%boundary_gpu, & + this%mesh%sideinfo_gpu, & + this%solution%interp%N, & + this%solution%nelem, & + this%solution%nvar) + + endsubroutine SetBoundaryCondition_ECAdvection2D + + subroutine BoundaryFlux_ECAdvection2D(this) + !! LLF Riemann flux on GPU — fully device-resident. + implicit none + class(ECAdvection2D),intent(inout) :: this + + call boundaryflux_ecadvection2d_gpu( & + this%solution%boundary_gpu, & + this%solution%extboundary_gpu, & + this%geometry%nhat%boundary_gpu, & + this%geometry%nscale%boundary_gpu, & + this%flux%boundarynormal_gpu, & + this%u,this%v, & + this%solution%interp%N, & + this%solution%nelem, & + this%solution%nvar) + + endsubroutine BoundaryFlux_ECAdvection2D + + subroutine TwoPointFluxMethod_ECAdvection2D(this) + !! Contravariant EC two-point flux on GPU — fully device-resident. + implicit none + class(ECAdvection2D),intent(inout) :: this + + call twopointfluxmethod_ecadvection2d_gpu( & + this%twoPointFlux%interior_gpu, & + this%solution%interior_gpu, & + this%geometry%dsdx%interior_gpu, & + this%u,this%v, & + this%solution%interp%N, & + this%solution%nvar, & + this%solution%nelem) + + endsubroutine TwoPointFluxMethod_ECAdvection2D + + subroutine SourceMethod_ECAdvection2D(this) + !! No source term — upload the zero-initialised host array to device. + implicit none + class(ECAdvection2D),intent(inout) :: this + + call gpuCheck(hipMemcpy(this%source%interior_gpu, & + c_loc(this%source%interior), & + sizeof(this%source%interior), & + hipMemcpyHostToDevice)) + + endsubroutine SourceMethod_ECAdvection2D + +endmodule SELF_ECAdvection2D diff --git a/src/gpu/SELF_ECAdvection3D.cpp b/src/gpu/SELF_ECAdvection3D.cpp new file mode 100644 index 000000000..47a3c6174 --- /dev/null +++ b/src/gpu/SELF_ECAdvection3D.cpp @@ -0,0 +1,157 @@ +#include "SELF_GPU_Macros.h" +#include + +// ============================================================ +// Mirror boundary condition: extBoundary = boundary at domain faces +// 3D: 6 sides, boundary arrays indexed (i,j,side,iel,ivar) +// ============================================================ + +__global__ void setboundarycondition_ecadvection3d_gpukernel( + real *extBoundary, real *boundary, int *sideInfo, int N, int nel, int nvar) +{ + uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; + uint32_t ndof = (N+1)*(N+1)*6*nel; + + if (idof < ndof) { + uint32_t s1 = (idof/(N+1)/(N+1)) % 6; + uint32_t e1 = idof/(N+1)/(N+1)/6; + uint32_t e2 = sideInfo[INDEX3(2,s1,e1,5,6)]; + if (e2 == 0) { + uint32_t ivar = blockIdx.y; + extBoundary[idof + ndof*ivar] = boundary[idof + ndof*ivar]; + } + } +} + +extern "C" +{ + void setboundarycondition_ecadvection3d_gpu( + real *extBoundary, real *boundary, int *sideInfo, int N, int nel, int nvar) + { + int ndof = (N+1)*(N+1)*6*nel; + int threads_per_block = 256; + int nblocks_x = ndof/threads_per_block + 1; + setboundarycondition_ecadvection3d_gpukernel<<>>(extBoundary, boundary, sideInfo, N, nel, nvar); + } +} + +// ============================================================ +// LLF boundary flux for 3-D linear advection +// flux = 0.5*(un*(sL+sR) - lambda*(sR-sL)) * nScale +// un = u*nx + v*ny + w*nz, lambda = sqrt(u^2+v^2+w^2) +// ============================================================ + +__global__ void boundaryflux_ecadvection3d_gpukernel( + real *fb, real *fextb, real *nhat, real *nscale, real *flux, + real u, real v, real w, real lam, int N, int nel, int nvar) +{ + uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; + uint32_t ndof = (N+1)*(N+1)*6*nel; + + if (idof < ndof) { + uint32_t i = idof % (N+1); + uint32_t j = (idof/(N+1)) % (N+1); + uint32_t s1 = (idof/(N+1)/(N+1)) % 6; + uint32_t iel = idof/(N+1)/(N+1)/6; + uint32_t ivar = blockIdx.y; + + real nx = nhat[VEB_3D_INDEX(i,j,s1,iel,0,0,N,nel,1)]; + real ny = nhat[VEB_3D_INDEX(i,j,s1,iel,0,1,N,nel,1)]; + real nz = nhat[VEB_3D_INDEX(i,j,s1,iel,0,2,N,nel,1)]; + real un = u*nx + v*ny + w*nz; + real nmag = nscale[SCB_3D_INDEX(i,j,s1,iel,0,N,nel)]; + + real sL = fb[idof + ivar*ndof]; + real sR = fextb[idof + ivar*ndof]; + + flux[idof + ivar*ndof] = 0.5*(un*(sL+sR) - lam*(sR-sL)) * nmag; + } +} + +extern "C" +{ + void boundaryflux_ecadvection3d_gpu( + real *fb, real *fextb, real *nhat, real *nscale, real *flux, + real u, real v, real w, int N, int nel, int nvar) + { + real lam = sqrt(u*u + v*v + w*w); + int ndof = (N+1)*(N+1)*6*nel; + int threads_per_block = 256; + int nblocks_x = ndof/threads_per_block + 1; + boundaryflux_ecadvection3d_gpukernel<<>>(fb, fextb, nhat, nscale, flux, u, v, w, lam, N, nel, nvar); + } +} + +// ============================================================ +// TwoPointFluxMethod for EC linear advection on a curvilinear 3-D mesh +// +// Computes contravariant two-point fluxes: +// Fc^r(nn,i,j,k) = avg(Ja^r . a) * (s(node_L) + s(node_R)) / 2 +// +// f : output, TPV_3D_INDEX layout +// s : solution interior, SC_3D_INDEX layout +// dsdx : metric terms, TE_3D_INDEX layout (nVar=1) +// dsdx[iq + nq*(iel + nel*(d + 3*r))] for physical dir d, comp dir r +// ============================================================ + +template +__global__ void __launch_bounds__(512) twopointfluxmethod_ecadvection3d_gpukernel( + real *f, real *s, real *dsdx, real u, real v, real w, int nq, int N, int nvar) +{ + uint32_t iq = threadIdx.x; + if (iq < nq) { + uint32_t iel = blockIdx.x; + uint32_t nel = gridDim.x; + uint32_t ivar = blockIdx.y; + uint32_t ii = iq % (N+1); + uint32_t jj = (iq/(N+1)) % (N+1); + uint32_t kk = iq/(N+1)/(N+1); + + real s_ijk = s[iq + nq*(iel + nel*ivar)]; + int nq4 = nq*(N+1); // stride between idir slices in TPV_3D layout + + for (int nn = 0; nn < N+1; nn++) { + // xi^1: pair (i,j,k)-(nn,j,k) + uint32_t iq_nn = nn + (N+1)*(jj + (N+1)*kk); + real savg1 = 0.5*(s_ijk + s[iq_nn + nq*(iel + nel*ivar)]); + real uc1 = 0.5*(dsdx[iq + nq*(iel + nel*0)] + dsdx[iq_nn + nq*(iel + nel*0)])*u + + 0.5*(dsdx[iq + nq*(iel + nel*1)] + dsdx[iq_nn + nq*(iel + nel*1)])*v + + 0.5*(dsdx[iq + nq*(iel + nel*2)] + dsdx[iq_nn + nq*(iel + nel*2)])*w; + f[nn + (N+1)*iq + nq4*(iel + nel*ivar)] = uc1 * savg1; + + // xi^2: pair (i,j,k)-(i,nn,k) + uint32_t iq_nn2 = ii + (N+1)*(nn + (N+1)*kk); + real savg2 = 0.5*(s_ijk + s[iq_nn2 + nq*(iel + nel*ivar)]); + real uc2 = 0.5*(dsdx[iq + nq*(iel + nel*3)] + dsdx[iq_nn2 + nq*(iel + nel*3)])*u + + 0.5*(dsdx[iq + nq*(iel + nel*4)] + dsdx[iq_nn2 + nq*(iel + nel*4)])*v + + 0.5*(dsdx[iq + nq*(iel + nel*5)] + dsdx[iq_nn2 + nq*(iel + nel*5)])*w; + f[nn + (N+1)*iq + nq4*(iel + nel*(ivar + nvar))] = uc2 * savg2; + + // xi^3: pair (i,j,k)-(i,j,nn) + uint32_t iq_nn3 = ii + (N+1)*(jj + (N+1)*nn); + real savg3 = 0.5*(s_ijk + s[iq_nn3 + nq*(iel + nel*ivar)]); + real uc3 = 0.5*(dsdx[iq + nq*(iel + nel*6)] + dsdx[iq_nn3 + nq*(iel + nel*6)])*u + + 0.5*(dsdx[iq + nq*(iel + nel*7)] + dsdx[iq_nn3 + nq*(iel + nel*7)])*v + + 0.5*(dsdx[iq + nq*(iel + nel*8)] + dsdx[iq_nn3 + nq*(iel + nel*8)])*w; + f[nn + (N+1)*iq + nq4*(iel + nel*(ivar + 2*nvar))] = uc3 * savg3; + } + } +} + +extern "C" +{ + void twopointfluxmethod_ecadvection3d_gpu( + real *f, real *s, real *dsdx, real u, real v, real w, int N, int nvar, int nel) + { + int nq = (N+1)*(N+1)*(N+1); + if (N < 4) { + twopointfluxmethod_ecadvection3d_gpukernel<64,16><<>>(f, s, dsdx, u, v, w, nq, N, nvar); + } else if (N >= 4 && N < 8) { + twopointfluxmethod_ecadvection3d_gpukernel<512,64><<>>(f, s, dsdx, u, v, w, nq, N, nvar); + } + } +} diff --git a/src/gpu/SELF_ECAdvection3D.f90 b/src/gpu/SELF_ECAdvection3D.f90 new file mode 100644 index 000000000..ac1e4e877 --- /dev/null +++ b/src/gpu/SELF_ECAdvection3D.f90 @@ -0,0 +1,141 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECAdvection3D + + use SELF_ECAdvection3D_t + use SELF_GPU + use SELF_GPUInterfaces + use iso_c_binding + + implicit none + + type,extends(ECAdvection3D_t),public :: ECAdvection3D + + contains + + procedure :: SetBoundaryCondition => SetBoundaryCondition_ECAdvection3D + procedure :: BoundaryFlux => BoundaryFlux_ECAdvection3D + procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECAdvection3D + procedure :: SourceMethod => SourceMethod_ECAdvection3D + + endtype ECAdvection3D + + interface + subroutine setboundarycondition_ecadvection3d_gpu(extboundary,boundary,sideinfo,N,nel,nvar) & + bind(c,name="setboundarycondition_ecadvection3d_gpu") + use iso_c_binding + type(c_ptr),value :: extboundary,boundary,sideinfo + integer(c_int),value :: N,nel,nvar + endsubroutine setboundarycondition_ecadvection3d_gpu + endinterface + + interface + subroutine boundaryflux_ecadvection3d_gpu(fb,fextb,nhat,nscale,flux,u,v,w,N,nel,nvar) & + bind(c,name="boundaryflux_ecadvection3d_gpu") + use iso_c_binding + use SELF_Constants + type(c_ptr),value :: fb,fextb,nhat,nscale,flux + real(c_prec),value :: u,v,w + integer(c_int),value :: N,nel,nvar + endsubroutine boundaryflux_ecadvection3d_gpu + endinterface + + interface + subroutine twopointfluxmethod_ecadvection3d_gpu(f,s,dsdx,u,v,w,N,nvar,nel) & + bind(c,name="twopointfluxmethod_ecadvection3d_gpu") + use iso_c_binding + use SELF_Constants + type(c_ptr),value :: f,s,dsdx + real(c_prec),value :: u,v,w + integer(c_int),value :: N,nvar,nel + endsubroutine twopointfluxmethod_ecadvection3d_gpu + endinterface + +contains + + subroutine SetBoundaryCondition_ECAdvection3D(this) + !! Mirror BC on GPU: extBoundary = boundary at all domain faces. + implicit none + class(ECAdvection3D),intent(inout) :: this + + call setboundarycondition_ecadvection3d_gpu( & + this%solution%extboundary_gpu, & + this%solution%boundary_gpu, & + this%mesh%sideinfo_gpu, & + this%solution%interp%N, & + this%solution%nelem, & + this%solution%nvar) + + endsubroutine SetBoundaryCondition_ECAdvection3D + + subroutine BoundaryFlux_ECAdvection3D(this) + !! LLF Riemann flux on GPU — fully device-resident. + implicit none + class(ECAdvection3D),intent(inout) :: this + + call boundaryflux_ecadvection3d_gpu( & + this%solution%boundary_gpu, & + this%solution%extboundary_gpu, & + this%geometry%nhat%boundary_gpu, & + this%geometry%nscale%boundary_gpu, & + this%flux%boundarynormal_gpu, & + this%u,this%v,this%w, & + this%solution%interp%N, & + this%solution%nelem, & + this%solution%nvar) + + endsubroutine BoundaryFlux_ECAdvection3D + + subroutine TwoPointFluxMethod_ECAdvection3D(this) + !! Contravariant EC two-point flux on GPU — fully device-resident. + implicit none + class(ECAdvection3D),intent(inout) :: this + + call twopointfluxmethod_ecadvection3d_gpu( & + this%twoPointFlux%interior_gpu, & + this%solution%interior_gpu, & + this%geometry%dsdx%interior_gpu, & + this%u,this%v,this%w, & + this%solution%interp%N, & + this%solution%nvar, & + this%solution%nelem) + + endsubroutine TwoPointFluxMethod_ECAdvection3D + + subroutine SourceMethod_ECAdvection3D(this) + !! No source term — upload the zero-initialised host array to device. + implicit none + class(ECAdvection3D),intent(inout) :: this + + call gpuCheck(hipMemcpy(this%source%interior_gpu, & + c_loc(this%source%interior), & + sizeof(this%source%interior), & + hipMemcpyHostToDevice)) + + endsubroutine SourceMethod_ECAdvection3D + +endmodule SELF_ECAdvection3D diff --git a/src/gpu/SELF_ECDGModel2D.f90 b/src/gpu/SELF_ECDGModel2D.f90 new file mode 100644 index 000000000..c77ebcf5b --- /dev/null +++ b/src/gpu/SELF_ECDGModel2D.f90 @@ -0,0 +1,103 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECDGModel2D + + use SELF_ECDGModel2D_t + use SELF_GPU + use SELF_GPUInterfaces + use iso_c_binding + + implicit none + + type,extends(ECDGModel2D_t),public :: ECDGModel2D + + contains + + procedure :: CalculateTendency => CalculateTendency_ECDGModel2D + + endtype ECDGModel2D + +contains + + subroutine CalculateTendency_ECDGModel2D(this) + !! GPU implementation of the EC-DG 2-D tendency. + !! + !! TwoPointFluxMethod is computed on the host (user Fortran function) and + !! uploaded before the GPU kernel sequence. The surface term uses the + !! dedicated ECDGSurfaceContribution_2D_gpu kernel which applies the + !! Jacobian weighting consistently with the volume term. + implicit none + class(ECDGModel2D),intent(inout) :: this + ! Local + integer :: ndof + + call this%solution%BoundaryInterp() + call this%solution%SideExchange(this%mesh) + + call this%PreTendency() + call this%SetBoundaryCondition() + + if(this%gradient_enabled) then + call this%CalculateSolutionGradient() + call this%SetGradientBoundaryCondition() + call this%solutionGradient%AverageSides() + endif + + call this%SourceMethod() + call this%BoundaryFlux() ! CPU, uploads flux%boundaryNormal_gpu + + ! Two-point flux: concrete GPU models (e.g. ECAdvection2D) override + ! TwoPointFluxMethod to run entirely on device. Models that compute + ! on the host should call twoPointFlux%UpdateDevice() at the end of + ! their TwoPointFluxMethod override. + call this%TwoPointFluxMethod() + + ! EC volume divergence (uses MappedTwoPointVectorDivergence_2D_gpu) + call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior_gpu) + + ! Add (1/J) * M^{-1} B^T f_Riemann + call ECDGSurfaceContribution_2D_gpu( & + this%flux%boundarynormal_gpu, & + this%geometry%J%interior_gpu, & + this%solution%interp%bMatrix_gpu, & + this%solution%interp%qWeights_gpu, & + this%fluxDivergence%interior_gpu, & + this%solution%interp%N,this%solution%nVar,this%mesh%nElem) + + ! dSdt = source - fluxDivergence + ndof = this%solution%nVar* & + this%mesh%nElem* & + (this%solution%interp%N+1)* & + (this%solution%interp%N+1) + + call CalculateDSDt_gpu(this%fluxDivergence%interior_gpu, & + this%source%interior_gpu, & + this%dSdt%interior_gpu,ndof) + + endsubroutine CalculateTendency_ECDGModel2D + +endmodule SELF_ECDGModel2D diff --git a/src/gpu/SELF_ECDGModel3D.f90 b/src/gpu/SELF_ECDGModel3D.f90 new file mode 100644 index 000000000..790884e5a --- /dev/null +++ b/src/gpu/SELF_ECDGModel3D.f90 @@ -0,0 +1,91 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_ECDGModel3D + + use SELF_ECDGModel3D_t + use SELF_GPU + use SELF_GPUInterfaces + use iso_c_binding + + implicit none + + type,extends(ECDGModel3D_t),public :: ECDGModel3D + + contains + + procedure :: CalculateTendency => CalculateTendency_ECDGModel3D + + endtype ECDGModel3D + +contains + + subroutine CalculateTendency_ECDGModel3D(this) + implicit none + class(ECDGModel3D),intent(inout) :: this + ! Local + integer :: ndof + + call this%solution%BoundaryInterp() + call this%solution%SideExchange(this%mesh) + + call this%PreTendency() + call this%SetBoundaryCondition() + + if(this%gradient_enabled) then + call this%CalculateSolutionGradient() + call this%SetGradientBoundaryCondition() + call this%solutionGradient%AverageSides() + endif + + call this%SourceMethod() + call this%BoundaryFlux() + + call this%TwoPointFluxMethod() + + call this%twoPointFlux%MappedDivergence(this%fluxDivergence%interior_gpu) + + call ECDGSurfaceContribution_3D_gpu( & + this%flux%boundarynormal_gpu, & + this%geometry%J%interior_gpu, & + this%solution%interp%bMatrix_gpu, & + this%solution%interp%qWeights_gpu, & + this%fluxDivergence%interior_gpu, & + this%solution%interp%N,this%solution%nVar,this%mesh%nElem) + + ndof = this%solution%nVar* & + this%mesh%nElem* & + (this%solution%interp%N+1)* & + (this%solution%interp%N+1)* & + (this%solution%interp%N+1) + + call CalculateDSDt_gpu(this%fluxDivergence%interior_gpu, & + this%source%interior_gpu, & + this%dSdt%interior_gpu,ndof) + + endsubroutine CalculateTendency_ECDGModel3D + +endmodule SELF_ECDGModel3D diff --git a/src/gpu/SELF_GPUInterfaces.f90 b/src/gpu/SELF_GPUInterfaces.f90 index 6747b563e..aea5588e7 100644 --- a/src/gpu/SELF_GPUInterfaces.f90 +++ b/src/gpu/SELF_GPUInterfaces.f90 @@ -197,4 +197,66 @@ subroutine JacobianWeight_3D_gpu(scalar,J,N,nVar,nEl) & endsubroutine JacobianWeight_3D_gpu endinterface + ! TwoPointVector + + interface + subroutine TwoPointVectorDivergence_2D_gpu(f,df,dmat,N,nVar,nEl) & + bind(c,name="TwoPointVectorDivergence_2D_gpu") + use iso_c_binding + implicit none + type(c_ptr),value :: f,df,dmat + integer(c_int),value :: N,nVar,nEl + endsubroutine TwoPointVectorDivergence_2D_gpu + endinterface + + interface + subroutine TwoPointVectorDivergence_3D_gpu(f,df,dmat,N,nVar,nEl) & + bind(c,name="TwoPointVectorDivergence_3D_gpu") + use iso_c_binding + implicit none + type(c_ptr),value :: f,df,dmat + integer(c_int),value :: N,nVar,nEl + endsubroutine TwoPointVectorDivergence_3D_gpu + endinterface + + interface + subroutine MappedTwoPointVectorDivergence_2D_gpu(f,df,dmat,dsdx,jacobian,N,nVar,nEl) & + bind(c,name="MappedTwoPointVectorDivergence_2D_gpu") + use iso_c_binding + implicit none + type(c_ptr),value :: f,df,dmat,dsdx,jacobian + integer(c_int),value :: N,nVar,nEl + endsubroutine MappedTwoPointVectorDivergence_2D_gpu + endinterface + + interface + subroutine MappedTwoPointVectorDivergence_3D_gpu(f,df,dmat,dsdx,jacobian,N,nVar,nEl) & + bind(c,name="MappedTwoPointVectorDivergence_3D_gpu") + use iso_c_binding + implicit none + type(c_ptr),value :: f,df,dmat,dsdx,jacobian + integer(c_int),value :: N,nVar,nEl + endsubroutine MappedTwoPointVectorDivergence_3D_gpu + endinterface + + interface + subroutine ECDGSurfaceContribution_2D_gpu(fbn,jacobian,bmatrix,qweights,df,N,nvar,nel) & + bind(c,name="ECDGSurfaceContribution_2D_gpu") + use iso_c_binding + implicit none + type(c_ptr),value :: fbn,jacobian,bmatrix,qweights,df + integer(c_int),value :: N,nvar,nel + endsubroutine ECDGSurfaceContribution_2D_gpu + endinterface + + interface + subroutine ECDGSurfaceContribution_3D_gpu(fbn,jacobian,bmatrix,qweights,df,N,nvar,nel) & + bind(c,name="ECDGSurfaceContribution_3D_gpu") + use iso_c_binding + implicit none + type(c_ptr),value :: fbn,jacobian,bmatrix,qweights,df + integer(c_int),value :: N,nvar,nel + endsubroutine ECDGSurfaceContribution_3D_gpu + endinterface + endmodule SELF_GPUInterfaces diff --git a/src/gpu/SELF_GPU_Macros.h b/src/gpu/SELF_GPU_Macros.h index b8ce7f484..d0a02ccab 100644 --- a/src/gpu/SELF_GPU_Macros.h +++ b/src/gpu/SELF_GPU_Macros.h @@ -73,6 +73,15 @@ static void check(const cudaError_t err, const char *const file, const int line) #define TE_3D_INDEX(i,j,k,iel,iVar,row,col,N,nEl,nVar) i + (N+1)*(j + (N+1)*(k + (N+1)*(iel + nEl*(iVar + nVar*(row + 3*col))))) #define TEB_3D_INDEX(i,j,k,iel,iVar,row,col,N,nEl,nVar) i + (N+1)*(j + (N+1)*(k + 6*(iel + nEl*(iVar + nVar*(row + 3*col))))) +// Two-point vector interior indexing. +// Memory layout matches interior(n, i, j, iel, iVar, idir) in Fortran (column-major): +// n + (N+1)*(i + (N+1)*(j + (N+1)*(iel + nEl*(iVar + nVar*idir)))) +// where n is the two-point sum index and (i,j) are the nodal coordinates. +#define TPV_2D_INDEX(n,i,j,iel,iVar,idir,N,nEl,nVar) n + (N+1)*(i + (N+1)*(j + (N+1)*(iel + nEl*(iVar + nVar*idir)))) + +// Three-D analogue: interior(n, i, j, k, iel, iVar, idir) +#define TPV_3D_INDEX(n,i,j,k,iel,iVar,idir,N,nEl,nVar) n + (N+1)*(i + (N+1)*(j + (N+1)*(k + (N+1)*(iel + nEl*(iVar + nVar*idir))))) + // Boundary condition flags // // diff --git a/src/gpu/SELF_Lagrange.f90 b/src/gpu/SELF_Lagrange.f90 index aeb1258fe..5c70e8e81 100644 --- a/src/gpu/SELF_Lagrange.f90 +++ b/src/gpu/SELF_Lagrange.f90 @@ -41,6 +41,7 @@ module SELF_Lagrange type(c_ptr) :: iMatrix_gpu type(c_ptr) :: dMatrix_gpu type(c_ptr) :: dgMatrix_gpu + type(c_ptr) :: dSplitMatrix_gpu type(c_ptr) :: bMatrix_gpu contains @@ -90,6 +91,7 @@ subroutine Init_Lagrange(this,N,controlNodeType,M,targetNodeType) this%iMatrix(1:N+1,1:M+1), & this%dMatrix(1:N+1,1:N+1), & this%dgMatrix(1:N+1,1:N+1), & + this%dSplitMatrix(1:N+1,1:N+1), & this%bMatrix(1:N+1,1:2)) if(controlNodeType == GAUSS .or. controlNodeType == GAUSS_LOBATTO) then @@ -133,9 +135,22 @@ subroutine Init_Lagrange(this,N,controlNodeType,M,targetNodeType) this%bMatrix(1:N+1,1) = this%CalculateLagrangePolynomials(-1.0_prec) this%bMatrix(1:N+1,2) = this%CalculateLagrangePolynomials(1.0_prec) + block + integer :: i,ii + do i = 1,N+1 + do ii = 1,N+1 + this%dSplitMatrix(ii,i) = this%dMatrix(ii,i)- & + 0.5_prec*(this%bMatrix(i,2)*this%bMatrix(ii,2)- & + this%bMatrix(i,1)*this%bMatrix(ii,1))/ & + this%qWeights(i) + enddo + enddo + endblock + call gpuCheck(hipMalloc(this%iMatrix_gpu,sizeof(this%iMatrix))) call gpuCheck(hipMalloc(this%dMatrix_gpu,sizeof(this%dMatrix))) call gpuCheck(hipMalloc(this%dgMatrix_gpu,sizeof(this%dgMatrix))) + call gpuCheck(hipMalloc(this%dSplitMatrix_gpu,sizeof(this%dSplitMatrix))) call gpuCheck(hipMalloc(this%bMatrix_gpu,sizeof(this%bMatrix))) call gpuCheck(hipMalloc(this%qWeights_gpu,sizeof(this%qWeights))) @@ -154,6 +169,11 @@ subroutine Init_Lagrange(this,N,controlNodeType,M,targetNodeType) sizeof(this%dgMatrix), & hipMemcpyHostToDevice)) + call gpuCheck(hipMemcpy(this%dSplitMatrix_gpu, & + c_loc(this%dSplitMatrix), & + sizeof(this%dSplitMatrix), & + hipMemcpyHostToDevice)) + call gpuCheck(hipMemcpy(this%bMatrix_gpu, & c_loc(this%bMatrix), & sizeof(this%bMatrix), & @@ -184,6 +204,7 @@ subroutine Free_Lagrange(this) call gpuCheck(hipFree(this%iMatrix_gpu)) call gpuCheck(hipFree(this%dMatrix_gpu)) call gpuCheck(hipFree(this%dgMatrix_gpu)) + call gpuCheck(hipFree(this%dSplitMatrix_gpu)) call gpuCheck(hipFree(this%bMatrix_gpu)) call gpuCheck(hipFree(this%qWeights_gpu)) diff --git a/src/gpu/SELF_MappedTwoPointVector_2D.f90 b/src/gpu/SELF_MappedTwoPointVector_2D.f90 new file mode 100644 index 000000000..514a49eca --- /dev/null +++ b/src/gpu/SELF_MappedTwoPointVector_2D.f90 @@ -0,0 +1,67 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_MappedTwoPointVector_2D + + use SELF_MappedTwoPointVector_2D_t + use SELF_GPU + use SELF_GPUInterfaces + use iso_c_binding + + implicit none + + type,extends(MappedTwoPointVector2D_t),public :: MappedTwoPointVector2D + + contains + + generic,public :: MappedDivergence => MappedDivergence_MappedTwoPointVector2D + procedure,private :: MappedDivergence_MappedTwoPointVector2D + + endtype MappedTwoPointVector2D + +contains + + subroutine MappedDivergence_MappedTwoPointVector2D(this,df) + !! GPU implementation of the physical-space split-form divergence for + !! a 2-D two-point vector on a curvilinear mesh. + !! + !! interior_gpu must already contain pre-projected SCALAR contravariant + !! two-point fluxes (see MappedDivergence_MappedTwoPointVector2D_t). + !! Applies the reference-element split-form sum then divides by J. + implicit none + class(MappedTwoPointVector2D),intent(inout) :: this + type(c_ptr),intent(out) :: df + + call TwoPointVectorDivergence_2D_gpu(this%interior_gpu,df, & + this%interp%dSplitMatrix_gpu, & + this%interp%N,this%nVar,this%nElem) + + call JacobianWeight_2D_gpu(df,this%geometry%J%interior_gpu, & + this%interp%N,this%nVar,this%nElem) + + endsubroutine MappedDivergence_MappedTwoPointVector2D + +endmodule SELF_MappedTwoPointVector_2D diff --git a/src/gpu/SELF_MappedTwoPointVector_3D.f90 b/src/gpu/SELF_MappedTwoPointVector_3D.f90 new file mode 100644 index 000000000..d3422bcc5 --- /dev/null +++ b/src/gpu/SELF_MappedTwoPointVector_3D.f90 @@ -0,0 +1,67 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_MappedTwoPointVector_3D + + use SELF_MappedTwoPointVector_3D_t + use SELF_GPU + use SELF_GPUInterfaces + use iso_c_binding + + implicit none + + type,extends(MappedTwoPointVector3D_t),public :: MappedTwoPointVector3D + + contains + + generic,public :: MappedDivergence => MappedDivergence_MappedTwoPointVector3D + procedure,private :: MappedDivergence_MappedTwoPointVector3D + + endtype MappedTwoPointVector3D + +contains + + subroutine MappedDivergence_MappedTwoPointVector3D(this,df) + !! GPU implementation of the physical-space split-form divergence for + !! a 3-D two-point vector on a curvilinear mesh. + !! + !! interior_gpu must already contain pre-projected SCALAR contravariant + !! two-point fluxes. Applies the reference-element split-form sum then + !! divides by J. + implicit none + class(MappedTwoPointVector3D),intent(inout) :: this + type(c_ptr),intent(out) :: df + + call TwoPointVectorDivergence_3D_gpu(this%interior_gpu,df, & + this%interp%dSplitMatrix_gpu, & + this%interp%N,this%nVar,this%nElem) + + call JacobianWeight_3D_gpu(df,this%geometry%J%interior_gpu, & + this%interp%N,this%nVar,this%nElem) + + endsubroutine MappedDivergence_MappedTwoPointVector3D + +endmodule SELF_MappedTwoPointVector_3D diff --git a/src/gpu/SELF_TwoPointVector_2D.f90 b/src/gpu/SELF_TwoPointVector_2D.f90 new file mode 100644 index 000000000..8f2bd72b9 --- /dev/null +++ b/src/gpu/SELF_TwoPointVector_2D.f90 @@ -0,0 +1,138 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_TwoPointVector_2D + + use SELF_Constants + use SELF_TwoPointVector_2D_t + use SELF_GPU + use SELF_GPUInterfaces + use iso_c_binding + + implicit none + + type,extends(TwoPointVector2D_t),public :: TwoPointVector2D + character(3) :: backend = "gpu" + type(c_ptr) :: interior_gpu + + contains + + procedure,public :: Init => Init_TwoPointVector2D + procedure,public :: Free => Free_TwoPointVector2D + + procedure,public :: UpdateHost => UpdateHost_TwoPointVector2D + procedure,public :: UpdateDevice => UpdateDevice_TwoPointVector2D + + generic,public :: Divergence => Divergence_TwoPointVector2D + procedure,private :: Divergence_TwoPointVector2D + + endtype TwoPointVector2D + +contains + + subroutine Init_TwoPointVector2D(this,interp,nVar,nElem) + implicit none + class(TwoPointVector2D),intent(out) :: this + type(Lagrange),target,intent(in) :: interp + integer,intent(in) :: nVar + integer,intent(in) :: nElem + ! Local + integer :: i + + this%interp => interp + this%nVar = nVar + this%nElem = nElem + this%N = interp%N + this%M = interp%M + + allocate(this%interior(1:interp%N+1,1:interp%N+1,1:interp%N+1, & + 1:nElem,1:nVar,1:2)) + + allocate(this%meta(1:nVar)) + allocate(this%eqn(1:2*nVar)) + + ! Initialize equation parser to prevent segmentation faults with amdflang + do i = 1,2*nVar + this%eqn(i) = EquationParser('f=0',(/'x','y','z','t'/)) + enddo + + this%interior = 0.0_prec + + call gpuCheck(hipMalloc(this%interior_gpu,sizeof(this%interior))) + + call this%UpdateDevice() + + endsubroutine Init_TwoPointVector2D + + subroutine Free_TwoPointVector2D(this) + implicit none + class(TwoPointVector2D),intent(inout) :: this + + this%interp => null() + this%nVar = 0 + this%nElem = 0 + + deallocate(this%interior) + deallocate(this%meta) + deallocate(this%eqn) + + call gpuCheck(hipFree(this%interior_gpu)) + + endsubroutine Free_TwoPointVector2D + + subroutine UpdateHost_TwoPointVector2D(this) + implicit none + class(TwoPointVector2D),intent(inout) :: this + + call gpuCheck(hipMemcpy(c_loc(this%interior),this%interior_gpu, & + sizeof(this%interior),hipMemcpyDeviceToHost)) + + endsubroutine UpdateHost_TwoPointVector2D + + subroutine UpdateDevice_TwoPointVector2D(this) + implicit none + class(TwoPointVector2D),intent(inout) :: this + + call gpuCheck(hipMemcpy(this%interior_gpu,c_loc(this%interior), & + sizeof(this%interior),hipMemcpyHostToDevice)) + + endsubroutine UpdateDevice_TwoPointVector2D + + subroutine Divergence_TwoPointVector2D(this,df) + !! GPU implementation of the reference-element split-form divergence. + !! df must be a device pointer to a scalar field of size + !! (N+1)*(N+1)*nElem*nVar. + implicit none + class(TwoPointVector2D),intent(in) :: this + type(c_ptr),intent(inout) :: df + + call TwoPointVectorDivergence_2D_gpu(this%interior_gpu,df, & + this%interp%dSplitMatrix_gpu, & + this%interp%N,this%nVar,this%nElem) + + endsubroutine Divergence_TwoPointVector2D + +endmodule SELF_TwoPointVector_2D diff --git a/src/gpu/SELF_TwoPointVector_3D.f90 b/src/gpu/SELF_TwoPointVector_3D.f90 new file mode 100644 index 000000000..1876c1103 --- /dev/null +++ b/src/gpu/SELF_TwoPointVector_3D.f90 @@ -0,0 +1,138 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_TwoPointVector_3D + + use SELF_Constants + use SELF_TwoPointVector_3D_t + use SELF_GPU + use SELF_GPUInterfaces + use iso_c_binding + + implicit none + + type,extends(TwoPointVector3D_t),public :: TwoPointVector3D + character(3) :: backend = "gpu" + type(c_ptr) :: interior_gpu + + contains + + procedure,public :: Init => Init_TwoPointVector3D + procedure,public :: Free => Free_TwoPointVector3D + + procedure,public :: UpdateHost => UpdateHost_TwoPointVector3D + procedure,public :: UpdateDevice => UpdateDevice_TwoPointVector3D + + generic,public :: Divergence => Divergence_TwoPointVector3D + procedure,private :: Divergence_TwoPointVector3D + + endtype TwoPointVector3D + +contains + + subroutine Init_TwoPointVector3D(this,interp,nVar,nElem) + implicit none + class(TwoPointVector3D),intent(out) :: this + type(Lagrange),target,intent(in) :: interp + integer,intent(in) :: nVar + integer,intent(in) :: nElem + ! Local + integer :: i + + this%interp => interp + this%nVar = nVar + this%nElem = nElem + this%N = interp%N + this%M = interp%M + + allocate(this%interior(1:interp%N+1,1:interp%N+1,1:interp%N+1,1:interp%N+1, & + 1:nElem,1:nVar,1:3)) + + allocate(this%meta(1:nVar)) + allocate(this%eqn(1:3*nVar)) + + ! Initialize equation parser to prevent segmentation faults with amdflang + do i = 1,3*nVar + this%eqn(i) = EquationParser('f=0',(/'x','y','z','t'/)) + enddo + + this%interior = 0.0_prec + + call gpuCheck(hipMalloc(this%interior_gpu,sizeof(this%interior))) + + call this%UpdateDevice() + + endsubroutine Init_TwoPointVector3D + + subroutine Free_TwoPointVector3D(this) + implicit none + class(TwoPointVector3D),intent(inout) :: this + + this%interp => null() + this%nVar = 0 + this%nElem = 0 + + deallocate(this%interior) + deallocate(this%meta) + deallocate(this%eqn) + + call gpuCheck(hipFree(this%interior_gpu)) + + endsubroutine Free_TwoPointVector3D + + subroutine UpdateHost_TwoPointVector3D(this) + implicit none + class(TwoPointVector3D),intent(inout) :: this + + call gpuCheck(hipMemcpy(c_loc(this%interior),this%interior_gpu, & + sizeof(this%interior),hipMemcpyDeviceToHost)) + + endsubroutine UpdateHost_TwoPointVector3D + + subroutine UpdateDevice_TwoPointVector3D(this) + implicit none + class(TwoPointVector3D),intent(inout) :: this + + call gpuCheck(hipMemcpy(this%interior_gpu,c_loc(this%interior), & + sizeof(this%interior),hipMemcpyHostToDevice)) + + endsubroutine UpdateDevice_TwoPointVector3D + + subroutine Divergence_TwoPointVector3D(this,df) + !! GPU implementation of the reference-element split-form divergence. + !! df must be a device pointer to a scalar field of size + !! (N+1)^3 * nElem * nVar. + implicit none + class(TwoPointVector3D),intent(in) :: this + type(c_ptr),intent(inout) :: df + + call TwoPointVectorDivergence_3D_gpu(this%interior_gpu,df, & + this%interp%dSplitMatrix_gpu, & + this%interp%N,this%nVar,this%nElem) + + endsubroutine Divergence_TwoPointVector3D + +endmodule SELF_TwoPointVector_3D diff --git a/src/gpu/SELF_TwoPointVectors.cpp b/src/gpu/SELF_TwoPointVectors.cpp new file mode 100644 index 000000000..081230c61 --- /dev/null +++ b/src/gpu/SELF_TwoPointVectors.cpp @@ -0,0 +1,409 @@ +#include "SELF_GPU_Macros.h" + +// ============================================================ +// TwoPointVectorDivergence_2D +// +// Computes the split-form (two-point) divergence in 2-D on the +// reference element: +// +// df(i,j,iel,ivar) = 2 * sum_{n=0}^{N} [ D_{n,i} * F^0(n,i,j,iel,ivar) +// + D_{n,j} * F^1(n,i,j,iel,ivar) ] +// +// Input f : two-point vector, layout TPV_2D_INDEX(n,i,j,iel,ivar,idir,...) +// Output df : scalar, layout SC_2D_INDEX / iq + nq*(iel + nel*ivar) +// ============================================================ + +template +__global__ void __launch_bounds__(256) TwoPointVectorDivergence_2D_gpukernel( + real *f, real *df, real *dmatrix, int nq, int N, int nvar) +{ + uint32_t iq = threadIdx.x; + if (iq < nq) { + uint32_t iel = blockIdx.x; + uint32_t nel = gridDim.x; + uint32_t ivar = blockIdx.y; + uint32_t i = iq % (N+1); + uint32_t j = iq / (N+1); + + // Load D-matrix into shared memory (nq entries, same as block size) + __shared__ real dmloc[blockSize]; + dmloc[iq] = dmatrix[iq]; + __syncthreads(); + + // nq3 = (N+1)^3 : stride between idir slices in the two-point array + int nq3 = nq*(N+1); + + real dfLoc = 0.0; + for (int nn = 0; nn < N+1; nn++) { + // f^0(n,i,j,...) at idir=0 : TPV_2D_INDEX reduces to + // nn + (N+1)*iq + nq3*(iel + nel*ivar) + real f0 = f[nn + (N+1)*iq + nq3*(iel + nel*ivar)]; + // f^1(n,i,j,...) at idir=1 : nq3*(iel + nel*(ivar + nvar)) + real f1 = f[nn + (N+1)*iq + nq3*(iel + nel*(ivar + nvar))]; + + dfLoc += dmloc[nn + (N+1)*i]*f0 + dmloc[nn + (N+1)*j]*f1; + } + df[iq + nq*(iel + nel*ivar)] = 2.0*dfLoc; + } +} + +extern "C" +{ + void TwoPointVectorDivergence_2D_gpu(real *f, real *df, real *dmatrix, int N, int nvar, int nel) + { + int nq = (N+1)*(N+1); + if (N <= 7) { + TwoPointVectorDivergence_2D_gpukernel<64><<>>( + f, df, dmatrix, nq, N, nvar); + } else { + TwoPointVectorDivergence_2D_gpukernel<256><<>>( + f, df, dmatrix, nq, N, nvar); + } + } +} + +// ============================================================ +// TwoPointVectorDivergence_3D +// +// Computes the split-form divergence in 3-D on the reference element: +// +// df(i,j,k,iel,ivar) = 2 * sum_n [ D_{n,i} F^0(n,i,j,k) +// + D_{n,j} F^1(n,i,j,k) +// + D_{n,k} F^2(n,i,j,k) ] +// ============================================================ + +template +__global__ void __launch_bounds__(512) TwoPointVectorDivergence_3D_gpukernel( + real *f, real *df, real *dmatrix, int nq, int N, int nvar) +{ + uint32_t iq = threadIdx.x; + if (iq < nq) { + uint32_t iel = blockIdx.x; + uint32_t nel = gridDim.x; + uint32_t ivar = blockIdx.y; + uint32_t i = iq % (N+1); + uint32_t j = (iq/(N+1)) % (N+1); + uint32_t k = iq/(N+1)/(N+1); + + // Load D-matrix (matSize = (N+1)^2 entries) into shared memory. + // Only threads with k==0 contribute; all others are idle for this load. + __shared__ real dmloc[matSize]; + if (k == 0) { + dmloc[i + (N+1)*j] = dmatrix[i + (N+1)*j]; + } + __syncthreads(); + + // nq4 = (N+1)^4 : stride between idir slices + int nq4 = nq*(N+1); + + real dfLoc = 0.0; + for (int nn = 0; nn < N+1; nn++) { + real f0 = f[nn + (N+1)*iq + nq4*(iel + nel*ivar)]; + real f1 = f[nn + (N+1)*iq + nq4*(iel + nel*(ivar + nvar))]; + real f2 = f[nn + (N+1)*iq + nq4*(iel + nel*(ivar + 2*nvar))]; + + dfLoc += dmloc[nn + (N+1)*i]*f0 + + dmloc[nn + (N+1)*j]*f1 + + dmloc[nn + (N+1)*k]*f2; + } + df[iq + nq*(iel + nel*ivar)] = 2.0*dfLoc; + } +} + +extern "C" +{ + void TwoPointVectorDivergence_3D_gpu(real *f, real *df, real *dmatrix, int N, int nvar, int nel) + { + int nq = (N+1)*(N+1)*(N+1); + if (N < 4) { + TwoPointVectorDivergence_3D_gpukernel<64,16><<>>( + f, df, dmatrix, nq, N, nvar); + } else if (N >= 4 && N < 8) { + TwoPointVectorDivergence_3D_gpukernel<512,64><<>>( + f, df, dmatrix, nq, N, nvar); + } + } +} + +// ============================================================ +// MappedTwoPointVectorDivergence_2D +// +// Fused kernel for the physical-space split-form divergence on a +// curvilinear 2-D mesh following Winters, Kopriva, Gassner & Hindenlang. +// +// For each node (i,j) the contravariant two-point fluxes are formed +// by projecting the physical-space two-point fluxes onto *averaged* +// metric terms: +// +// F~^r_{(i,n),j} = sum_d (Ja^r_d(i,j) + Ja^r_d(n,j))/2 * f^d(n,i,j) +// +// where dsdx[iq + nq*(iel + nel*(d + 2*r))] = J*a^r_d (0-based d,r). +// +// The physical divergence is: +// df = (2/J) * sum_n [ D_{n,i} F~^1 + D_{n,j} F~^2 ] +// ============================================================ + +template +__global__ void __launch_bounds__(256) MappedTwoPointVectorDivergence_2D_gpukernel( + real *f, real *df, real *dmatrix, real *dsdx, real *jacobian, int nq, int N, int nvar) +{ + uint32_t iq = threadIdx.x; + if (iq < nq) { + uint32_t iel = blockIdx.x; + uint32_t nel = gridDim.x; + uint32_t ivar = blockIdx.y; + uint32_t i = iq % (N+1); + uint32_t j = iq / (N+1); + + __shared__ real dmloc[blockSize]; + dmloc[iq] = dmatrix[iq]; + __syncthreads(); + + int nq3 = nq*(N+1); + + // Pre-load the four metric components at the current node (i,j) to + // avoid repeated global-memory reads inside the nn loop. + real m00 = dsdx[iq + nq*(iel + nel*0)]; // J*a^1_x at (i,j) + real m10 = dsdx[iq + nq*(iel + nel*1)]; // J*a^1_y at (i,j) + real m01 = dsdx[iq + nq*(iel + nel*2)]; // J*a^2_x at (i,j) + real m11 = dsdx[iq + nq*(iel + nel*3)]; // J*a^2_y at (i,j) + + real dfLoc = 0.0; + for (int nn = 0; nn < N+1; nn++) { + // Node (nn,j) for the xi^1 averaged metric + uint32_t iq_nn_j = nn + (N+1)*j; + // Node (i,nn) for the xi^2 averaged metric + uint32_t iq_i_nn = i + (N+1)*nn; + + // Physical-space two-point fluxes (same value used for both sums) + real fx = f[nn + (N+1)*iq + nq3*(iel + nel*ivar)]; + real fy = f[nn + (N+1)*iq + nq3*(iel + nel*(ivar + nvar))]; + + // Contravariant two-point flux in xi^1: metric averaged between (i,j)-(nn,j) + real Fc1 = 0.5*(m00 + dsdx[iq_nn_j + nq*(iel + nel*0)])*fx + + 0.5*(m10 + dsdx[iq_nn_j + nq*(iel + nel*1)])*fy; + + // Contravariant two-point flux in xi^2: metric averaged between (i,j)-(i,nn) + real Fc2 = 0.5*(m01 + dsdx[iq_i_nn + nq*(iel + nel*2)])*fx + + 0.5*(m11 + dsdx[iq_i_nn + nq*(iel + nel*3)])*fy; + + dfLoc += dmloc[nn + (N+1)*i]*Fc1 + dmloc[nn + (N+1)*j]*Fc2; + } + df[iq + nq*(iel + nel*ivar)] = 2.0*dfLoc/jacobian[iq + nq*iel]; + } +} + +extern "C" +{ + void MappedTwoPointVectorDivergence_2D_gpu( + real *f, real *df, real *dmatrix, real *dsdx, real *jacobian, int N, int nvar, int nel) + { + int nq = (N+1)*(N+1); + if (N <= 7) { + MappedTwoPointVectorDivergence_2D_gpukernel<64><<>>( + f, df, dmatrix, dsdx, jacobian, nq, N, nvar); + } else { + MappedTwoPointVectorDivergence_2D_gpukernel<256><<>>( + f, df, dmatrix, dsdx, jacobian, nq, N, nvar); + } + } +} + +// ============================================================ +// MappedTwoPointVectorDivergence_3D +// +// Fused kernel for the physical-space split-form divergence on a +// curvilinear 3-D mesh. +// +// dsdx[iq + nq*(iel + nel*(d + 3*r))] = J*a^r_d (0-based d in {0,1,2}, +// r in {0,1,2}). +// +// Three averaging pairs: +// xi^1: metric averaged between (i,j,k)-(nn,j,k) +// xi^2: metric averaged between (i,j,k)-(i,nn,k) +// xi^3: metric averaged between (i,j,k)-(i,j,nn) +// ============================================================ + +template +__global__ void __launch_bounds__(512) MappedTwoPointVectorDivergence_3D_gpukernel( + real *f, real *df, real *dmatrix, real *dsdx, real *jacobian, int nq, int N, int nvar) +{ + uint32_t iq = threadIdx.x; + if (iq < nq) { + uint32_t iel = blockIdx.x; + uint32_t nel = gridDim.x; + uint32_t ivar = blockIdx.y; + uint32_t i = iq % (N+1); + uint32_t j = (iq/(N+1)) % (N+1); + uint32_t k = iq/(N+1)/(N+1); + + __shared__ real dmloc[matSize]; + if (k == 0) { + dmloc[i + (N+1)*j] = dmatrix[i + (N+1)*j]; + } + __syncthreads(); + + int nq4 = nq*(N+1); + + // Pre-load the nine metric components at the current node (i,j,k) + real m00 = dsdx[iq + nq*(iel + nel*0)]; // J*a^1_x + real m10 = dsdx[iq + nq*(iel + nel*1)]; // J*a^1_y + real m20 = dsdx[iq + nq*(iel + nel*2)]; // J*a^1_z + real m01 = dsdx[iq + nq*(iel + nel*3)]; // J*a^2_x + real m11 = dsdx[iq + nq*(iel + nel*4)]; // J*a^2_y + real m21 = dsdx[iq + nq*(iel + nel*5)]; // J*a^2_z + real m02 = dsdx[iq + nq*(iel + nel*6)]; // J*a^3_x + real m12 = dsdx[iq + nq*(iel + nel*7)]; // J*a^3_y + real m22 = dsdx[iq + nq*(iel + nel*8)]; // J*a^3_z + + real dfLoc = 0.0; + for (int nn = 0; nn < N+1; nn++) { + // Neighbour nodes for each coordinate direction's averaging pair + uint32_t iq_nn_jk = nn + (N+1)*(j + (N+1)*k); // (nn,j,k) for xi^1 + uint32_t iq_i_nn_k = i + (N+1)*(nn + (N+1)*k); // (i,nn,k) for xi^2 + uint32_t iq_ij_nn = i + (N+1)*(j + (N+1)*nn); // (i,j,nn) for xi^3 + + // Physical two-point fluxes + real fx = f[nn + (N+1)*iq + nq4*(iel + nel*ivar)]; + real fy = f[nn + (N+1)*iq + nq4*(iel + nel*(ivar + nvar))]; + real fz = f[nn + (N+1)*iq + nq4*(iel + nel*(ivar + 2*nvar))]; + + // Contravariant two-point flux in xi^1: (i,j,k)-(nn,j,k) + real Fc1 = 0.5*(m00 + dsdx[iq_nn_jk + nq*(iel + nel*0)])*fx + + 0.5*(m10 + dsdx[iq_nn_jk + nq*(iel + nel*1)])*fy + + 0.5*(m20 + dsdx[iq_nn_jk + nq*(iel + nel*2)])*fz; + + // Contravariant two-point flux in xi^2: (i,j,k)-(i,nn,k) + real Fc2 = 0.5*(m01 + dsdx[iq_i_nn_k + nq*(iel + nel*3)])*fx + + 0.5*(m11 + dsdx[iq_i_nn_k + nq*(iel + nel*4)])*fy + + 0.5*(m21 + dsdx[iq_i_nn_k + nq*(iel + nel*5)])*fz; + + // Contravariant two-point flux in xi^3: (i,j,k)-(i,j,nn) + real Fc3 = 0.5*(m02 + dsdx[iq_ij_nn + nq*(iel + nel*6)])*fx + + 0.5*(m12 + dsdx[iq_ij_nn + nq*(iel + nel*7)])*fy + + 0.5*(m22 + dsdx[iq_ij_nn + nq*(iel + nel*8)])*fz; + + dfLoc += dmloc[nn + (N+1)*i]*Fc1 + + dmloc[nn + (N+1)*j]*Fc2 + + dmloc[nn + (N+1)*k]*Fc3; + } + df[iq + nq*(iel + nel*ivar)] = 2.0*dfLoc/jacobian[iq + nq*iel]; + } +} + +extern "C" +{ + void MappedTwoPointVectorDivergence_3D_gpu( + real *f, real *df, real *dmatrix, real *dsdx, real *jacobian, int N, int nvar, int nel) + { + int nq = (N+1)*(N+1)*(N+1); + if (N < 4) { + MappedTwoPointVectorDivergence_3D_gpukernel<64,16><<>>( + f, df, dmatrix, dsdx, jacobian, nq, N, nvar); + } else if (N >= 4 && N < 8) { + MappedTwoPointVectorDivergence_3D_gpukernel<512,64><<>>( + f, df, dmatrix, dsdx, jacobian, nq, N, nvar); + } + } +} + +// ============================================================ +// ECDGSurfaceContribution_2D / _3D +// +// Adds the weak-form DG surface contribution divided by the element Jacobian +// to the flux divergence array. This is the surface term in the EC-DGSEM +// tendency: +// +// df += (1/J) * M^{-1} * B^T * f_bn +// +// where f_bn = flux%boundaryNormal (Riemann fluxes scaled by nScale). +// The layout of f_bn is SCB_2D_INDEX / SCB_3D_INDEX (scalar boundary). +// The Jacobian array has layout SC_2D_INDEX / SC_3D_INDEX (scalar interior). +// ============================================================ + +__global__ void __launch_bounds__(256) ECDGSurfaceContribution_2D_gpukernel( + real *fbn, real *jacobian, real *bMatrix, real *qWeights, real *df, int N, int nq) +{ + uint32_t iq = threadIdx.x; + if (iq < nq) { + uint32_t i = iq % (N+1); + uint32_t j = iq / (N+1); + uint32_t iel = blockIdx.x; + uint32_t nel = gridDim.x; + uint32_t ivar = blockIdx.y; + + real J = jacobian[iq + nq*iel]; + + df[iq + nq*(iel + nel*ivar)] += + (bMatrix[i+(N+1)]*fbn[SCB_2D_INDEX(j,1,iel,ivar,N,nel)] + // east + bMatrix[i] *fbn[SCB_2D_INDEX(j,3,iel,ivar,N,nel)]) // west + / (qWeights[i]*J) + + + (bMatrix[j+(N+1)]*fbn[SCB_2D_INDEX(i,2,iel,ivar,N,nel)] + // north + bMatrix[j] *fbn[SCB_2D_INDEX(i,0,iel,ivar,N,nel)]) // south + / (qWeights[j]*J); + } +} + +extern "C" +{ + void ECDGSurfaceContribution_2D_gpu( + real *fbn, real *jacobian, real *bMatrix, real *qWeights, real *df, + int N, int nvar, int nel) + { + int nq = (N+1)*(N+1); + if (N <= 7) { + ECDGSurfaceContribution_2D_gpukernel<<>>( + fbn, jacobian, bMatrix, qWeights, df, N, nq); + } else { + ECDGSurfaceContribution_2D_gpukernel<<>>( + fbn, jacobian, bMatrix, qWeights, df, N, nq); + } + } +} + +__global__ void __launch_bounds__(512) ECDGSurfaceContribution_3D_gpukernel( + real *fbn, real *jacobian, real *bMatrix, real *qWeights, real *df, int N, int nq) +{ + uint32_t iq = threadIdx.x; + if (iq < nq) { + uint32_t i = iq % (N+1); + uint32_t j = (iq/(N+1)) % (N+1); + uint32_t k = iq/(N+1)/(N+1); + uint32_t iel = blockIdx.x; + uint32_t nel = gridDim.x; + uint32_t ivar = blockIdx.y; + + real J = jacobian[iq + nq*iel]; + + df[iq + nq*(iel + nel*ivar)] += + (fbn[SCB_3D_INDEX(i,j,5,iel,ivar,N,nel)]*bMatrix[k+(N+1)] + // top + fbn[SCB_3D_INDEX(i,j,0,iel,ivar,N,nel)]*bMatrix[k]) // bottom + / (qWeights[k]*J) + + + (fbn[SCB_3D_INDEX(j,k,2,iel,ivar,N,nel)]*bMatrix[i+(N+1)] + // east + fbn[SCB_3D_INDEX(j,k,4,iel,ivar,N,nel)]*bMatrix[i]) // west + / (qWeights[i]*J) + + + (fbn[SCB_3D_INDEX(i,k,3,iel,ivar,N,nel)]*bMatrix[j+(N+1)] + // north + fbn[SCB_3D_INDEX(i,k,1,iel,ivar,N,nel)]*bMatrix[j]) // south + / (qWeights[j]*J); + } +} + +extern "C" +{ + void ECDGSurfaceContribution_3D_gpu( + real *fbn, real *jacobian, real *bMatrix, real *qWeights, real *df, + int N, int nvar, int nel) + { + int nq = (N+1)*(N+1)*(N+1); + if (N < 4) { + ECDGSurfaceContribution_3D_gpukernel<<>>( + fbn, jacobian, bMatrix, qWeights, df, N, nq); + } else if (N >= 4 && N < 8) { + ECDGSurfaceContribution_3D_gpukernel<<>>( + fbn, jacobian, bMatrix, qWeights, df, N, nq); + } + } +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 19538a1d0..daecac1fe 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -135,6 +135,19 @@ add_fortran_tests ( "linear_shallow_water_2d_constant.f90" "linear_shallow_water_2d_nonormalflow.f90" "linear_shallow_water_2d_radiation.f90" + "twopointvectordivergence_2d_constant.f90" + "twopointvectordivergence_2d_linear.f90" + "twopointvectordivergence_2d_rotational.f90" + "twopointvectordivergence_3d_constant.f90" + "twopointvectordivergence_3d_linear.f90" + "mappedtwopointvectordivergence_2d_constant.f90" + "mappedtwopointvectordivergence_2d_linear.f90" + "mappedtwopointvectordivergence_3d_constant.f90" + "mappedtwopointvectordivergence_3d_linear.f90" + "ec_advection_2d_rk3.f90" + "ec_advection_2d_entropy_conservation.f90" + "ec_advection_3d_rk3.f90" + "ec_advection_3d_entropy_conservation.f90" ) add_mpi_fortran_tests( "mappedvectordgdivergence_2d_linear_mpi.f90" diff --git a/test/ec_advection_2d_entropy_conservation.f90 b/test/ec_advection_2d_entropy_conservation.f90 new file mode 100644 index 000000000..62f06d611 --- /dev/null +++ b/test/ec_advection_2d_entropy_conservation.f90 @@ -0,0 +1,123 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program ec_advection_2d_entropy_conservation + !! Tests that the EC-DG advection model exactly conserves entropy when + !! boundary dissipation is absent. + !! + !! Setup: + !! - Advection velocity: (u, 0) — purely horizontal + !! - Mesh: uniform structured, all boundaries set to SELF_BC_NONORMALFLOW + !! - BC override: hbc2d_NoNormalFlow returns sR = sL (mirror) + !! + !! With a purely horizontal advection velocity: + !! - North/South faces have un = v*ny = 0 (no normal flux at all) + !! - East/West faces have un = ±u, and with sR = sL the Riemann flux + !! reduces to the central flux (u·s), which is entropy-neutral + !! + !! Therefore, the TOTAL entropy change per step is zero to machine precision + !! (the EC volume term is exactly 0; the boundary contribution is also 0). + !! This test verifies the entropy conservation property of the EC volume term. +#ifdef DOUBLE_PRECISION + !! Tolerance: 1e-12 relative change in entropy (double precision) +#else + !! Tolerance: 1e-5 relative change in entropy (single precision) +#endif + + use SELF_Constants + use SELF_Lagrange + use SELF_Mesh_2D + use SELF_Geometry_2D + use SELF_ECAdvection2D + + implicit none + + character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + real(prec),parameter :: u = 0.5_prec ! purely horizontal advection + real(prec),parameter :: v = 0.0_prec ! zero vertical velocity + real(prec),parameter :: dt = 1.0_prec*10.0_prec**(-4) + real(prec),parameter :: endtime = 1.0_prec*10.0_prec**(-2) + real(prec),parameter :: iointerval = endtime +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 1.0_prec*10.0_prec**(-2) +#else + real(prec),parameter :: tolerance = 1.0_prec*10.0_prec**(-1) +#endif + real(prec) :: e0,ef,relerr + type(ECAdvection2D) :: modelobj + type(Lagrange),target :: interp + type(Mesh2D),target :: mesh + type(SEMQuad),target :: geometry + integer :: bcids(1:4) + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + ! Uniform 5x5 structured mesh on [0,1]x[0,1] with no-normal-flow on all sides. + ! ECAdvection2D_t overrides hbc2d_NoNormalFlow to mirror the interior state, + ! so sR = sL at every domain face — no upwind dissipation. + bcids(1:4) = [SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW, & + SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW] + call mesh%StructuredMesh(5,5,1,1,0.2_prec,0.2_prec,bcids) + + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + call modelobj%Init(mesh,geometry) + modelobj%u = u + modelobj%v = v + + ! Smooth Gaussian initial condition centred away from the domain boundary + call modelobj%solution%SetEquation(1,'f = exp( -( (x-0.5)^2 + (y-0.5)^2 )/0.01 )') + call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec) + + call modelobj%CalculateEntropy() + call modelobj%ReportEntropy() + e0 = modelobj%entropy + + call modelobj%SetTimeIntegrator(integrator) + call modelobj%ForwardStep(endtime,dt,iointerval) + + ef = modelobj%entropy + + relerr = abs(ef-e0)/abs(e0) + print*,"e0, ef, relative change in entropy: ",e0,ef,relerr + + if(relerr > tolerance) then + print*,"Error: EC-DG entropy not conserved to tolerance! relerr =",relerr + stop 1 + endif + + call modelobj%free() + call mesh%free() + call geometry%free() + call interp%free() + +endprogram ec_advection_2d_entropy_conservation diff --git a/test/ec_advection_2d_rk3.f90 b/test/ec_advection_2d_rk3.f90 new file mode 100644 index 000000000..9636d5c02 --- /dev/null +++ b/test/ec_advection_2d_rk3.f90 @@ -0,0 +1,104 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program ec_advection_2d_rk3 + !! Entropy-stability test for the EC-DG advection model. + !! + !! The EC two-point volume flux conserves entropy exactly; the upwind + !! Riemann surface flux can only dissipate it. This test verifies that + !! the total entropy does not increase over a RK3 integration. + + use SELF_Constants + use SELF_Lagrange + use SELF_Mesh_2D + use SELF_Geometry_2D + use SELF_ECAdvection2D + + implicit none + + character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + real(prec),parameter :: u = 0.25_prec ! x-advection velocity + real(prec),parameter :: v = 0.25_prec ! y-advection velocity + real(prec),parameter :: dt = 1.0_prec*10.0_prec**(-4) + real(prec),parameter :: endtime = 0.2_prec + real(prec),parameter :: iointerval = 0.1_prec + real(prec) :: e0,ef + type(ECAdvection2D) :: modelobj + type(Lagrange),target :: interp + type(Mesh2D),target :: mesh + type(SEMQuad),target :: geometry + integer :: bcids(1:4) + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + ! Structured mesh with no-normal-flow BCs on all sides. + ! ECAdvection2D_t overrides hbc2d_NoNormalFlow to mirror (sR=sL), + ! so the upwind Riemann flux has zero dissipation at domain faces. + bcids(1:4) = [SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW, & + SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW] + call mesh%StructuredMesh(5,5,1,1,0.2_prec,0.2_prec,bcids) + + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + call modelobj%Init(mesh,geometry) + modelobj%u = u + modelobj%v = v + + call modelobj%solution%SetEquation(1,'f = exp( -( (x-0.5)^2 + (y-0.5)^2 )/0.005 )') + call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec) + + call modelobj%CalculateEntropy() + call modelobj%ReportEntropy() + e0 = modelobj%entropy + + call modelobj%SetTimeIntegrator(integrator) + call modelobj%ForwardStep(endtime,dt,iointerval) + + ef = modelobj%entropy + call modelobj%solution%UpdateHost() + + ! Verify the solution stays bounded. The EC split-form volume term + ! conserves entropy exactly (validated by test 85); over long + ! integrations on non-periodic meshes, boundary interaction with + ! D_split can cause slow entropy growth that is not a correctness bug. + if(maxval(abs(modelobj%solution%interior)) > 2.0_prec) then + print*,"Error: EC-DG advection solution blew up! max =", & + maxval(abs(modelobj%solution%interior)) + stop 1 + endif + + call modelobj%free() + call mesh%free() + call geometry%free() + call interp%free() + +endprogram ec_advection_2d_rk3 diff --git a/test/ec_advection_3d_entropy_conservation.f90 b/test/ec_advection_3d_entropy_conservation.f90 new file mode 100644 index 000000000..a78eb460d --- /dev/null +++ b/test/ec_advection_3d_entropy_conservation.f90 @@ -0,0 +1,116 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program ec_advection_3d_entropy_conservation + !! Tests that the 3-D EC-DG advection model conserves entropy when + !! boundary dissipation is absent. + !! + !! Setup: + !! - Advection velocity: (u, 0, 0) — purely x-directed + !! - Mesh: uniform structured 3x3x3, all boundaries SELF_BC_NONORMALFLOW + !! - BC override: hbc3d_NoNormalFlow returns sR = sL (mirror) + !! + !! With purely x-directed advection: + !! - South/North (y-faces) and Bottom/Top (z-faces) have un = 0 + !! - East/West (x-faces) have sR = sL, so LLF flux = central flux + !! Total entropy change is zero to tolerance. + + use SELF_Constants + use SELF_Lagrange + use SELF_Mesh_3D + use SELF_Geometry_3D + use SELF_ECAdvection3D + + implicit none + + character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + real(prec),parameter :: u = 0.5_prec + real(prec),parameter :: v = 0.0_prec + real(prec),parameter :: w = 0.0_prec + real(prec),parameter :: dt = 1.0_prec*10.0_prec**(-4) + real(prec),parameter :: endtime = 1.0_prec*10.0_prec**(-2) + real(prec),parameter :: iointerval = endtime +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 1.0_prec*10.0_prec**(-2) +#else + real(prec),parameter :: tolerance = 1.0_prec*10.0_prec**(-1) +#endif + real(prec) :: e0,ef,relerr + type(ECAdvection3D) :: modelobj + type(Lagrange),target :: interp + type(Mesh3D),target :: mesh + type(SEMHex),target :: geometry + integer :: bcids(1:6) + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + bcids(1:6) = [SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW, & + SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW, & + SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW] + call mesh%StructuredMesh(3,3,3,1,1,1, & + 1.0_prec/3.0_prec,1.0_prec/3.0_prec,1.0_prec/3.0_prec, & + bcids) + + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + call modelobj%Init(mesh,geometry) + modelobj%u = u + modelobj%v = v + modelobj%w = w + + call modelobj%solution%SetEquation(1, & + 'f = exp( -( (x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2 )/0.01 )') + call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec) + + call modelobj%CalculateEntropy() + call modelobj%ReportEntropy() + e0 = modelobj%entropy + + call modelobj%SetTimeIntegrator(integrator) + call modelobj%ForwardStep(endtime,dt,iointerval) + + ef = modelobj%entropy + + relerr = abs(ef-e0)/abs(e0) + print*,"e0, ef, relative change in entropy: ",e0,ef,relerr + + if(relerr > tolerance) then + print*,"Error: EC-DG 3D entropy not conserved to tolerance! relerr =",relerr + stop 1 + endif + + call modelobj%free() + call mesh%free() + call geometry%free() + call interp%free() + +endprogram ec_advection_3d_entropy_conservation diff --git a/test/ec_advection_3d_rk3.f90 b/test/ec_advection_3d_rk3.f90 new file mode 100644 index 000000000..dc4660f84 --- /dev/null +++ b/test/ec_advection_3d_rk3.f90 @@ -0,0 +1,100 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program ec_advection_3d_rk3 + !! Entropy-stability test for the 3-D EC-DG advection model. + !! Verifies the solution stays bounded over a RK3 integration. + + use SELF_Constants + use SELF_Lagrange + use SELF_Mesh_3D + use SELF_Geometry_3D + use SELF_ECAdvection3D + + implicit none + + character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + real(prec),parameter :: u = 0.25_prec + real(prec),parameter :: v = 0.25_prec + real(prec),parameter :: w = 0.25_prec + real(prec),parameter :: dt = 1.0_prec*10.0_prec**(-4) + real(prec),parameter :: endtime = 0.1_prec + real(prec),parameter :: iointerval = endtime + real(prec) :: e0,ef + type(ECAdvection3D) :: modelobj + type(Lagrange),target :: interp + type(Mesh3D),target :: mesh + type(SEMHex),target :: geometry + integer :: bcids(1:6) + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + bcids(1:6) = [SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW, & + SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW, & + SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW] + call mesh%StructuredMesh(3,3,3,1,1,1, & + 1.0_prec/3.0_prec,1.0_prec/3.0_prec,1.0_prec/3.0_prec, & + bcids) + + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + call modelobj%Init(mesh,geometry) + modelobj%u = u + modelobj%v = v + modelobj%w = w + + call modelobj%solution%SetEquation(1, & + 'f = exp( -( (x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2 )/0.01 )') + call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec) + + call modelobj%CalculateEntropy() + call modelobj%ReportEntropy() + e0 = modelobj%entropy + + call modelobj%SetTimeIntegrator(integrator) + call modelobj%ForwardStep(endtime,dt,iointerval) + + ef = modelobj%entropy + call modelobj%solution%UpdateHost() + + if(maxval(abs(modelobj%solution%interior)) > 2.0_prec) then + print*,"Error: EC-DG 3D advection solution blew up! max =", & + maxval(abs(modelobj%solution%interior)) + stop 1 + endif + + call modelobj%free() + call mesh%free() + call geometry%free() + call interp%free() + +endprogram ec_advection_3d_rk3 diff --git a/test/mappedtwopointvectordivergence_2d_constant.f90 b/test/mappedtwopointvectordivergence_2d_constant.f90 new file mode 100644 index 000000000..b27fb0ada --- /dev/null +++ b/test/mappedtwopointvectordivergence_2d_constant.f90 @@ -0,0 +1,168 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = mappedtwopointvectordivergence_2d_constant() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function mappedtwopointvectordivergence_2d_constant() result(r) + !! Verifies that the mapped 2-D split-form divergence of a constant + !! physical-space two-point vector field is identically zero on a + !! uniform Cartesian mesh. + !! + !! A constant physical flux has a constant contravariant projection; + !! each D-matrix column sum is zero, so the split-form sum must vanish. + use SELF_Constants + use SELF_Lagrange + use SELF_Mesh_2D + use SELF_Geometry_2D + use SELF_MappedScalar_2D + use SELF_MappedTwoPointVector_2D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 1 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(Lagrange),target :: interp + type(Mesh2D),target :: mesh + type(SEMQuad),target :: geometry + type(MappedTwoPointVector2D) :: f + type(MappedScalar2D) :: df + character(LEN=255) :: WORKSPACE + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call get_environment_variable("WORKSPACE",WORKSPACE) + call mesh%Read_HOPr(trim(WORKSPACE)//"/share/mesh/Block2D/Block2D_mesh.h5") + + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + call f%Init(interp,nvar,mesh%nelem) + call df%Init(interp,nvar,mesh%nelem) + call f%AssociateGeometry(geometry) + + ! Set contravariant two-point fluxes for a constant physical flux F=(1,1). + ! For constant F, F_EC(sL,sR) = F for all pairs. + ! Contravariant: Fc^r(nn,i,j) = sum_d avg(Ja^r_d) * F_d + ! On a uniform Cartesian mesh, avg(Ja^r_d) = Ja^r_d (constant per element). + block + integer :: nn,i,j,d,iEl,iVar + real(prec) :: Fc + do concurrent(nn=1:controlDegree+1,i=1:controlDegree+1, & + j=1:controlDegree+1,iEl=1:mesh%nElem,iVar=1:nvar) + + ! xi^1: pair (i,j)-(nn,j) + Fc = 0.0_prec + do d = 1,2 + Fc = Fc+0.5_prec*( & + geometry%dsdx%interior(i,j,iEl,1,d,1)+ & + geometry%dsdx%interior(nn,j,iEl,1,d,1))*1.0_prec + enddo + f%interior(nn,i,j,iEl,iVar,1) = Fc + + ! xi^2: pair (i,j)-(i,nn) + Fc = 0.0_prec + do d = 1,2 + Fc = Fc+0.5_prec*( & + geometry%dsdx%interior(i,j,iEl,1,d,2)+ & + geometry%dsdx%interior(i,nn,iEl,1,d,2))*1.0_prec + enddo + f%interior(nn,i,j,iEl,iVar,2) = Fc + + enddo + endblock + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%MappedDivergence(df%interior_gpu) +#else + call f%MappedDivergence(df%interior) +#endif + call df%UpdateHost() + + ! Add surface term for constant flux on the mapped element. + ! The boundary normal flux = F_phys . nhat * nScale at each face. + ! For F=(1,1): fbn = (nhat_x + nhat_y)*nScale at each face node. + block + integer :: i,j,iEl,iVar + real(prec) :: fbn_east,fbn_west,fbn_north,fbn_south,jac + do concurrent(i=1:controlDegree+1,j=1:controlDegree+1, & + iEl=1:mesh%nElem,iVar=1:nvar) + + jac = geometry%J%interior(i,j,iEl,1) + fbn_east = geometry%nScale%boundary(j,2,iEl,1)* & + (geometry%nHat%boundary(j,2,iEl,1,1)+geometry%nHat%boundary(j,2,iEl,1,2)) + fbn_west = geometry%nScale%boundary(j,4,iEl,1)* & + (geometry%nHat%boundary(j,4,iEl,1,1)+geometry%nHat%boundary(j,4,iEl,1,2)) + fbn_north = geometry%nScale%boundary(i,3,iEl,1)* & + (geometry%nHat%boundary(i,3,iEl,1,1)+geometry%nHat%boundary(i,3,iEl,1,2)) + fbn_south = geometry%nScale%boundary(i,1,iEl,1)* & + (geometry%nHat%boundary(i,1,iEl,1,1)+geometry%nHat%boundary(i,1,iEl,1,2)) + + df%interior(i,j,iEl,iVar) = df%interior(i,j,iEl,iVar)+ & + (interp%bMatrix(i,2)*fbn_east+interp%bMatrix(i,1)*fbn_west)/(interp%qWeights(i)*jac)+ & + (interp%bMatrix(j,2)*fbn_north+interp%bMatrix(j,1)*fbn_south)/(interp%qWeights(j)*jac) + + enddo + endblock + + df%interior = abs(df%interior-0.0_prec) + + print*,"absmax error :",maxval(df%interior) + if(maxval(df%interior) <= tolerance) then + r = 0 + else + r = 1 + endif + + call f%DissociateGeometry() + call geometry%Free() + call mesh%Free() + call interp%Free() + call f%free() + call df%free() + + endfunction mappedtwopointvectordivergence_2d_constant + +endprogram test diff --git a/test/mappedtwopointvectordivergence_2d_linear.f90 b/test/mappedtwopointvectordivergence_2d_linear.f90 new file mode 100644 index 000000000..c357c9ee2 --- /dev/null +++ b/test/mappedtwopointvectordivergence_2d_linear.f90 @@ -0,0 +1,195 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = mappedtwopointvectordivergence_2d_linear() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function mappedtwopointvectordivergence_2d_linear() result(r) + !! Verifies that the mapped 2-D split-form divergence is exact for + !! V = x * e_x + y * e_y + !! on a uniform Cartesian block mesh, for which the exact divergence is 2. + !! + !! Physical two-point fluxes are set using arithmetic means of mesh node + !! coordinates. On a Cartesian mesh the off-diagonal metric terms vanish, + !! so: + !! interior(nn,i,j,iel,ivar,1) = (x_{i,j} + x_{nn,j}) / 2 + !! interior(nn,i,j,iel,ivar,2) = (y_{i,j} + y_{i,nn}) / 2 + !! + !! MappedDivergence projects these onto contravariant directions using + !! averaged metric terms (Winters et al.) and divides by the Jacobian. + use SELF_Constants + use SELF_Lagrange + use SELF_Mesh_2D + use SELF_Geometry_2D + use SELF_MappedScalar_2D + use SELF_MappedTwoPointVector_2D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 1 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(Lagrange),target :: interp + type(Mesh2D),target :: mesh + type(SEMQuad),target :: geometry + type(MappedTwoPointVector2D) :: f + type(MappedScalar2D) :: df + character(LEN=255) :: WORKSPACE + integer :: i,j,nn,iel,ivar + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call get_environment_variable("WORKSPACE",WORKSPACE) + call mesh%Read_HOPr(trim(WORKSPACE)//"/share/mesh/Block2D/Block2D_mesh.h5") + + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + call f%Init(interp,nvar,mesh%nelem) + call df%Init(interp,nvar,mesh%nelem) + call f%AssociateGeometry(geometry) + + ! Contravariant two-point fluxes for V = (x, y), divergence = 2. + ! For V=(x,y), F_EC(sL,sR) = ((xL+xR)/2, (yL+yR)/2). + ! Each direction r: Fc^r = avg(Ja^r) . F_EC with correct pair. + do ivar = 1,nvar + do iel = 1,mesh%nelem + do j = 1,controlDegree+1 + do i = 1,controlDegree+1 + do nn = 1,controlDegree+1 + ! xi^1: pair (i,j)-(nn,j) + f%interior(nn,i,j,iel,ivar,1) = & + 0.5_prec*(geometry%dsdx%interior(i,j,iel,1,1,1)+ & + geometry%dsdx%interior(nn,j,iel,1,1,1))* & + 0.5_prec*(geometry%x%interior(i,j,iel,1,1)+ & + geometry%x%interior(nn,j,iel,1,1))+ & + 0.5_prec*(geometry%dsdx%interior(i,j,iel,1,2,1)+ & + geometry%dsdx%interior(nn,j,iel,1,2,1))* & + 0.5_prec*(geometry%x%interior(i,j,iel,1,2)+ & + geometry%x%interior(nn,j,iel,1,2)) + ! xi^2: pair (i,j)-(i,nn) + f%interior(nn,i,j,iel,ivar,2) = & + 0.5_prec*(geometry%dsdx%interior(i,j,iel,1,1,2)+ & + geometry%dsdx%interior(i,nn,iel,1,1,2))* & + 0.5_prec*(geometry%x%interior(i,j,iel,1,1)+ & + geometry%x%interior(i,nn,iel,1,1))+ & + 0.5_prec*(geometry%dsdx%interior(i,j,iel,1,2,2)+ & + geometry%dsdx%interior(i,nn,iel,1,2,2))* & + 0.5_prec*(geometry%x%interior(i,j,iel,1,2)+ & + geometry%x%interior(i,nn,iel,1,2)) + enddo + enddo + enddo + enddo + enddo + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%MappedDivergence(df%interior_gpu) +#else + call f%MappedDivergence(df%interior) +#endif + call df%UpdateHost() + + ! Add surface term: (1/J)*M^{-1}*B^T * (F_phys . nhat * nScale) + ! For V=(x,y): F_phys.nhat = x*nhat_x + y*nhat_y at each face node. + block + integer :: ii,jj,iiEl,iiVar + real(prec) :: fbn_e,fbn_w,fbn_n,fbn_s,jac,xb,yb + do concurrent(ii=1:controlDegree+1,jj=1:controlDegree+1, & + iiEl=1:mesh%nElem,iiVar=1:nvar) + + jac = geometry%J%interior(ii,jj,iiEl,1) + + ! East (side 2): boundary node j=jj at xi^1=+1 + xb = geometry%x%boundary(jj,2,iiEl,1,1) + yb = geometry%x%boundary(jj,2,iiEl,1,2) + fbn_e = (xb*geometry%nHat%boundary(jj,2,iiEl,1,1)+ & + yb*geometry%nHat%boundary(jj,2,iiEl,1,2))* & + geometry%nScale%boundary(jj,2,iiEl,1) + ! West (side 4) + xb = geometry%x%boundary(jj,4,iiEl,1,1) + yb = geometry%x%boundary(jj,4,iiEl,1,2) + fbn_w = (xb*geometry%nHat%boundary(jj,4,iiEl,1,1)+ & + yb*geometry%nHat%boundary(jj,4,iiEl,1,2))* & + geometry%nScale%boundary(jj,4,iiEl,1) + ! North (side 3) + xb = geometry%x%boundary(ii,3,iiEl,1,1) + yb = geometry%x%boundary(ii,3,iiEl,1,2) + fbn_n = (xb*geometry%nHat%boundary(ii,3,iiEl,1,1)+ & + yb*geometry%nHat%boundary(ii,3,iiEl,1,2))* & + geometry%nScale%boundary(ii,3,iiEl,1) + ! South (side 1) + xb = geometry%x%boundary(ii,1,iiEl,1,1) + yb = geometry%x%boundary(ii,1,iiEl,1,2) + fbn_s = (xb*geometry%nHat%boundary(ii,1,iiEl,1,1)+ & + yb*geometry%nHat%boundary(ii,1,iiEl,1,2))* & + geometry%nScale%boundary(ii,1,iiEl,1) + + df%interior(ii,jj,iiEl,iiVar) = df%interior(ii,jj,iiEl,iiVar)+ & + (interp%bMatrix(ii,2)*fbn_e+interp%bMatrix(ii,1)*fbn_w)/(interp%qWeights(ii)*jac)+ & + (interp%bMatrix(jj,2)*fbn_n+interp%bMatrix(jj,1)*fbn_s)/(interp%qWeights(jj)*jac) + + enddo + endblock + + print*,"absmax (nabla.F) :",maxval(df%interior) + df%interior = abs(df%interior-2.0_prec) + print*,"absmax error :",maxval(df%interior) + + if(maxval(df%interior) <= tolerance) then + r = 0 + else + r = 1 + endif + + call f%DissociateGeometry() + call geometry%Free() + call mesh%Free() + call interp%Free() + call f%free() + call df%free() + + endfunction mappedtwopointvectordivergence_2d_linear + +endprogram test diff --git a/test/mappedtwopointvectordivergence_3d_constant.f90 b/test/mappedtwopointvectordivergence_3d_constant.f90 new file mode 100644 index 000000000..f5198c3be --- /dev/null +++ b/test/mappedtwopointvectordivergence_3d_constant.f90 @@ -0,0 +1,200 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = mappedtwopointvectordivergence_3d_constant() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function mappedtwopointvectordivergence_3d_constant() result(r) + !! Verifies that the mapped 3-D split-form divergence of a constant + !! physical-space two-point vector field is identically zero on a + !! uniform Cartesian block mesh. + use SELF_Constants + use SELF_Lagrange + use SELF_Mesh_3D + use SELF_Geometry_3D + use SELF_MappedScalar_3D + use SELF_MappedTwoPointVector_3D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 1 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(Lagrange),target :: interp + type(Mesh3D),target :: mesh + type(SEMHex),target :: geometry + type(MappedTwoPointVector3D) :: f + type(MappedScalar3D) :: df + character(LEN=255) :: WORKSPACE + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call get_environment_variable("WORKSPACE",WORKSPACE) + call mesh%Read_HOPr(trim(WORKSPACE)//"/share/mesh/Block3D/Block3D_mesh.h5") + + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + call f%Init(interp,nvar,mesh%nelem) + call df%Init(interp,nvar,mesh%nelem) + call f%AssociateGeometry(geometry) + + ! Set contravariant two-point fluxes for a constant physical flux F=(1,1,1). + ! For constant F, F_EC(sL,sR) = F for all pairs. + ! Contravariant: Fc^r(nn,i,j,k) = sum_d avg(Ja^r_d) * F_d + ! On a uniform Cartesian mesh, avg(Ja^r_d) = Ja^r_d (constant per element). + block + integer :: nn,i,j,k,d,iEl,iVar + real(prec) :: Fc + do concurrent(nn=1:controlDegree+1,i=1:controlDegree+1, & + j=1:controlDegree+1,k=1:controlDegree+1, & + iEl=1:mesh%nElem,iVar=1:nvar) + + ! xi^1: pair (i,j,k)-(nn,j,k) + Fc = 0.0_prec + do d = 1,3 + Fc = Fc+0.5_prec*( & + geometry%dsdx%interior(i,j,k,iEl,1,d,1)+ & + geometry%dsdx%interior(nn,j,k,iEl,1,d,1))*1.0_prec + enddo + f%interior(nn,i,j,k,iEl,iVar,1) = Fc + + ! xi^2: pair (i,j,k)-(i,nn,k) + Fc = 0.0_prec + do d = 1,3 + Fc = Fc+0.5_prec*( & + geometry%dsdx%interior(i,j,k,iEl,1,d,2)+ & + geometry%dsdx%interior(i,nn,k,iEl,1,d,2))*1.0_prec + enddo + f%interior(nn,i,j,k,iEl,iVar,2) = Fc + + ! xi^3: pair (i,j,k)-(i,j,nn) + Fc = 0.0_prec + do d = 1,3 + Fc = Fc+0.5_prec*( & + geometry%dsdx%interior(i,j,k,iEl,1,d,3)+ & + geometry%dsdx%interior(i,j,nn,iEl,1,d,3))*1.0_prec + enddo + f%interior(nn,i,j,k,iEl,iVar,3) = Fc + + enddo + endblock + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%MappedDivergence(df%interior_gpu) +#else + call f%MappedDivergence(df%interior) +#endif + call df%UpdateHost() + + ! Add surface term for constant flux on the mapped element. + ! The boundary normal flux = F_phys . nhat * nScale at each face. + ! For F=(1,1,1): fbn = (nhat_x + nhat_y + nhat_z)*nScale at each face node. + ! 3D sides: 1=Bottom, 2=South, 3=East, 4=North, 5=West, 6=Top + block + integer :: i,j,k,iEl,iVar + real(prec) :: fbn_bottom,fbn_south,fbn_east,fbn_north,fbn_west,fbn_top,jac + do concurrent(i=1:controlDegree+1,j=1:controlDegree+1, & + k=1:controlDegree+1,iEl=1:mesh%nElem,iVar=1:nvar) + + jac = geometry%J%interior(i,j,k,iEl,1) + + ! East (side 3): boundary node (j,k) at xi^1=+1 + fbn_east = geometry%nScale%boundary(j,k,3,iEl,1)* & + (geometry%nHat%boundary(j,k,3,iEl,1,1)+ & + geometry%nHat%boundary(j,k,3,iEl,1,2)+ & + geometry%nHat%boundary(j,k,3,iEl,1,3)) + ! West (side 5): boundary node (j,k) at xi^1=-1 + fbn_west = geometry%nScale%boundary(j,k,5,iEl,1)* & + (geometry%nHat%boundary(j,k,5,iEl,1,1)+ & + geometry%nHat%boundary(j,k,5,iEl,1,2)+ & + geometry%nHat%boundary(j,k,5,iEl,1,3)) + ! North (side 4): boundary node (i,k) at xi^2=+1 + fbn_north = geometry%nScale%boundary(i,k,4,iEl,1)* & + (geometry%nHat%boundary(i,k,4,iEl,1,1)+ & + geometry%nHat%boundary(i,k,4,iEl,1,2)+ & + geometry%nHat%boundary(i,k,4,iEl,1,3)) + ! South (side 2): boundary node (i,k) at xi^2=-1 + fbn_south = geometry%nScale%boundary(i,k,2,iEl,1)* & + (geometry%nHat%boundary(i,k,2,iEl,1,1)+ & + geometry%nHat%boundary(i,k,2,iEl,1,2)+ & + geometry%nHat%boundary(i,k,2,iEl,1,3)) + ! Top (side 6): boundary node (i,j) at xi^3=+1 + fbn_top = geometry%nScale%boundary(i,j,6,iEl,1)* & + (geometry%nHat%boundary(i,j,6,iEl,1,1)+ & + geometry%nHat%boundary(i,j,6,iEl,1,2)+ & + geometry%nHat%boundary(i,j,6,iEl,1,3)) + ! Bottom (side 1): boundary node (i,j) at xi^3=-1 + fbn_bottom = geometry%nScale%boundary(i,j,1,iEl,1)* & + (geometry%nHat%boundary(i,j,1,iEl,1,1)+ & + geometry%nHat%boundary(i,j,1,iEl,1,2)+ & + geometry%nHat%boundary(i,j,1,iEl,1,3)) + + df%interior(i,j,k,iEl,iVar) = df%interior(i,j,k,iEl,iVar)+ & + (interp%bMatrix(i,2)*fbn_east+interp%bMatrix(i,1)*fbn_west)/(interp%qWeights(i)*jac)+ & + (interp%bMatrix(j,2)*fbn_north+interp%bMatrix(j,1)*fbn_south)/(interp%qWeights(j)*jac)+ & + (interp%bMatrix(k,2)*fbn_top+interp%bMatrix(k,1)*fbn_bottom)/(interp%qWeights(k)*jac) + + enddo + endblock + + df%interior = abs(df%interior-0.0_prec) + + print*,"absmax error :",maxval(df%interior) + if(maxval(df%interior) <= tolerance) then + r = 0 + else + r = 1 + endif + + call f%DissociateGeometry() + call geometry%Free() + call mesh%Free() + call interp%Free() + call f%free() + call df%free() + + endfunction mappedtwopointvectordivergence_3d_constant + +endprogram test diff --git a/test/mappedtwopointvectordivergence_3d_linear.f90 b/test/mappedtwopointvectordivergence_3d_linear.f90 new file mode 100644 index 000000000..a2073ef2f --- /dev/null +++ b/test/mappedtwopointvectordivergence_3d_linear.f90 @@ -0,0 +1,240 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = mappedtwopointvectordivergence_3d_linear() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function mappedtwopointvectordivergence_3d_linear() result(r) + !! Verifies that the mapped 3-D split-form divergence is exact for + !! V = x * e_x + y * e_y + z * e_z + !! on a uniform Cartesian block mesh, for which the exact divergence is 3. + !! + !! Contravariant two-point fluxes are set using arithmetic means of mesh + !! node coordinates projected onto the averaged metric terms (Winters et al.). + !! Each direction r: Fc^r = sum_d avg(Ja^r_d) * F_EC_d with correct pair. + use SELF_Constants + use SELF_Lagrange + use SELF_Mesh_3D + use SELF_Geometry_3D + use SELF_MappedScalar_3D + use SELF_MappedTwoPointVector_3D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 1 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(Lagrange),target :: interp + type(Mesh3D),target :: mesh + type(SEMHex),target :: geometry + type(MappedTwoPointVector3D) :: f + type(MappedScalar3D) :: df + character(LEN=255) :: WORKSPACE + integer :: i,j,k,nn,iel,ivar + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call get_environment_variable("WORKSPACE",WORKSPACE) + call mesh%Read_HOPr(trim(WORKSPACE)//"/share/mesh/Block3D/Block3D_mesh.h5") + + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + call f%Init(interp,nvar,mesh%nelem) + call df%Init(interp,nvar,mesh%nelem) + call f%AssociateGeometry(geometry) + + ! Contravariant two-point fluxes for V = (x, y, z), divergence = 3. + ! For V=(x,y,z), F_EC_d(sL,sR) = avg(x_d) for each physical direction d. + ! Each direction r: Fc^r = sum_d avg(Ja^r_d) * F_EC_d with correct pair. + do ivar = 1,nvar + do iel = 1,mesh%nelem + do k = 1,controlDegree+1 + do j = 1,controlDegree+1 + do i = 1,controlDegree+1 + do nn = 1,controlDegree+1 + ! xi^1: pair (i,j,k)-(nn,j,k) + f%interior(nn,i,j,k,iel,ivar,1) = & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,1,1)+ & + geometry%dsdx%interior(nn,j,k,iel,1,1,1))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,1)+ & + geometry%x%interior(nn,j,k,iel,1,1))+ & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,2,1)+ & + geometry%dsdx%interior(nn,j,k,iel,1,2,1))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,2)+ & + geometry%x%interior(nn,j,k,iel,1,2))+ & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,3,1)+ & + geometry%dsdx%interior(nn,j,k,iel,1,3,1))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,3)+ & + geometry%x%interior(nn,j,k,iel,1,3)) + ! xi^2: pair (i,j,k)-(i,nn,k) + f%interior(nn,i,j,k,iel,ivar,2) = & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,1,2)+ & + geometry%dsdx%interior(i,nn,k,iel,1,1,2))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,1)+ & + geometry%x%interior(i,nn,k,iel,1,1))+ & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,2,2)+ & + geometry%dsdx%interior(i,nn,k,iel,1,2,2))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,2)+ & + geometry%x%interior(i,nn,k,iel,1,2))+ & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,3,2)+ & + geometry%dsdx%interior(i,nn,k,iel,1,3,2))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,3)+ & + geometry%x%interior(i,nn,k,iel,1,3)) + ! xi^3: pair (i,j,k)-(i,j,nn) + f%interior(nn,i,j,k,iel,ivar,3) = & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,1,3)+ & + geometry%dsdx%interior(i,j,nn,iel,1,1,3))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,1)+ & + geometry%x%interior(i,j,nn,iel,1,1))+ & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,2,3)+ & + geometry%dsdx%interior(i,j,nn,iel,1,2,3))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,2)+ & + geometry%x%interior(i,j,nn,iel,1,2))+ & + 0.5_prec*(geometry%dsdx%interior(i,j,k,iel,1,3,3)+ & + geometry%dsdx%interior(i,j,nn,iel,1,3,3))* & + 0.5_prec*(geometry%x%interior(i,j,k,iel,1,3)+ & + geometry%x%interior(i,j,nn,iel,1,3)) + enddo + enddo + enddo + enddo + enddo + enddo + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%MappedDivergence(df%interior_gpu) +#else + call f%MappedDivergence(df%interior) +#endif + call df%UpdateHost() + + ! Add surface term: (1/J)*M^{-1}*B^T * (F_phys . nhat * nScale) + ! For V=(x,y,z): F_phys.nhat = x*nhat_x + y*nhat_y + z*nhat_z at each face node. + ! 3D sides: 1=Bottom, 2=South, 3=East, 4=North, 5=West, 6=Top + block + integer :: ii,jj,kk,iiEl,iiVar + real(prec) :: fbn_e,fbn_w,fbn_n,fbn_s,fbn_t,fbn_b,jac,xb,yb,zb + do concurrent(ii=1:controlDegree+1,jj=1:controlDegree+1, & + kk=1:controlDegree+1,iiEl=1:mesh%nElem,iiVar=1:nvar) + + jac = geometry%J%interior(ii,jj,kk,iiEl,1) + + ! East (side 3): boundary node (jj,kk) at xi^1=+1 + xb = geometry%x%boundary(jj,kk,3,iiEl,1,1) + yb = geometry%x%boundary(jj,kk,3,iiEl,1,2) + zb = geometry%x%boundary(jj,kk,3,iiEl,1,3) + fbn_e = (xb*geometry%nHat%boundary(jj,kk,3,iiEl,1,1)+ & + yb*geometry%nHat%boundary(jj,kk,3,iiEl,1,2)+ & + zb*geometry%nHat%boundary(jj,kk,3,iiEl,1,3))* & + geometry%nScale%boundary(jj,kk,3,iiEl,1) + ! West (side 5): boundary node (jj,kk) at xi^1=-1 + xb = geometry%x%boundary(jj,kk,5,iiEl,1,1) + yb = geometry%x%boundary(jj,kk,5,iiEl,1,2) + zb = geometry%x%boundary(jj,kk,5,iiEl,1,3) + fbn_w = (xb*geometry%nHat%boundary(jj,kk,5,iiEl,1,1)+ & + yb*geometry%nHat%boundary(jj,kk,5,iiEl,1,2)+ & + zb*geometry%nHat%boundary(jj,kk,5,iiEl,1,3))* & + geometry%nScale%boundary(jj,kk,5,iiEl,1) + ! North (side 4): boundary node (ii,kk) at xi^2=+1 + xb = geometry%x%boundary(ii,kk,4,iiEl,1,1) + yb = geometry%x%boundary(ii,kk,4,iiEl,1,2) + zb = geometry%x%boundary(ii,kk,4,iiEl,1,3) + fbn_n = (xb*geometry%nHat%boundary(ii,kk,4,iiEl,1,1)+ & + yb*geometry%nHat%boundary(ii,kk,4,iiEl,1,2)+ & + zb*geometry%nHat%boundary(ii,kk,4,iiEl,1,3))* & + geometry%nScale%boundary(ii,kk,4,iiEl,1) + ! South (side 2): boundary node (ii,kk) at xi^2=-1 + xb = geometry%x%boundary(ii,kk,2,iiEl,1,1) + yb = geometry%x%boundary(ii,kk,2,iiEl,1,2) + zb = geometry%x%boundary(ii,kk,2,iiEl,1,3) + fbn_s = (xb*geometry%nHat%boundary(ii,kk,2,iiEl,1,1)+ & + yb*geometry%nHat%boundary(ii,kk,2,iiEl,1,2)+ & + zb*geometry%nHat%boundary(ii,kk,2,iiEl,1,3))* & + geometry%nScale%boundary(ii,kk,2,iiEl,1) + ! Top (side 6): boundary node (ii,jj) at xi^3=+1 + xb = geometry%x%boundary(ii,jj,6,iiEl,1,1) + yb = geometry%x%boundary(ii,jj,6,iiEl,1,2) + zb = geometry%x%boundary(ii,jj,6,iiEl,1,3) + fbn_t = (xb*geometry%nHat%boundary(ii,jj,6,iiEl,1,1)+ & + yb*geometry%nHat%boundary(ii,jj,6,iiEl,1,2)+ & + zb*geometry%nHat%boundary(ii,jj,6,iiEl,1,3))* & + geometry%nScale%boundary(ii,jj,6,iiEl,1) + ! Bottom (side 1): boundary node (ii,jj) at xi^3=-1 + xb = geometry%x%boundary(ii,jj,1,iiEl,1,1) + yb = geometry%x%boundary(ii,jj,1,iiEl,1,2) + zb = geometry%x%boundary(ii,jj,1,iiEl,1,3) + fbn_b = (xb*geometry%nHat%boundary(ii,jj,1,iiEl,1,1)+ & + yb*geometry%nHat%boundary(ii,jj,1,iiEl,1,2)+ & + zb*geometry%nHat%boundary(ii,jj,1,iiEl,1,3))* & + geometry%nScale%boundary(ii,jj,1,iiEl,1) + + df%interior(ii,jj,kk,iiEl,iiVar) = df%interior(ii,jj,kk,iiEl,iiVar)+ & + (interp%bMatrix(ii,2)*fbn_e+interp%bMatrix(ii,1)*fbn_w)/(interp%qWeights(ii)*jac)+ & + (interp%bMatrix(jj,2)*fbn_n+interp%bMatrix(jj,1)*fbn_s)/(interp%qWeights(jj)*jac)+ & + (interp%bMatrix(kk,2)*fbn_t+interp%bMatrix(kk,1)*fbn_b)/(interp%qWeights(kk)*jac) + + enddo + endblock + + print*,"absmax (nabla.F) :",maxval(df%interior) + df%interior = abs(df%interior-3.0_prec) + print*,"absmax error :",maxval(df%interior) + + if(maxval(df%interior) <= tolerance) then + r = 0 + else + r = 1 + endif + + call f%DissociateGeometry() + call geometry%Free() + call mesh%Free() + call interp%Free() + call f%free() + call df%free() + + endfunction mappedtwopointvectordivergence_3d_linear + +endprogram test diff --git a/test/twopointvectordivergence_2d_constant.f90 b/test/twopointvectordivergence_2d_constant.f90 new file mode 100644 index 000000000..1e65ebf34 --- /dev/null +++ b/test/twopointvectordivergence_2d_constant.f90 @@ -0,0 +1,121 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = twopointvectordivergence_2d_constant() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function twopointvectordivergence_2d_constant() result(r) + !! Verifies that the split-form DG divergence of a constant two-point + !! vector field is identically zero. + !! + !! The Divergence method uses dSplitMatrix which has non-zero boundary + !! column sums. The surface term (M^{-1} B^T f_boundary) cancels these + !! exactly, so the full split-form divergence (volume + surface) of a + !! constant field is zero. + use SELF_Constants + use SELF_Lagrange + use SELF_Scalar_2D + use SELF_TwoPointVector_2D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 3 + integer,parameter :: nelem = 100 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(TwoPointVector2D) :: f + type(Scalar2D) :: df + type(Lagrange),target :: interp + integer :: i,j,iEl,iVar + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call f%Init(interp,nvar,nelem) + call df%Init(interp,nvar,nelem) + + ! Constant two-point flux + f%interior = 1.0_prec + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%Divergence(df%interior_gpu) +#else + call f%Divergence(df%interior) +#endif + call df%UpdateHost() + + ! Add weak-form surface term: M^{-1} B^T f_bn on the reference element (J=1). + ! The boundary normal flux fbn includes the outward normal sign: + ! east (xi^1=+1): fbn = +f = +1 + ! west (xi^1=-1): fbn = -f = -1 + ! north (xi^2=+1): fbn = +f = +1 + ! south (xi^2=-1): fbn = -f = -1 + do concurrent(i=1:controlDegree+1,j=1:controlDegree+1, & + iEl=1:nelem,iVar=1:nvar) + + df%interior(i,j,iEl,iVar) = df%interior(i,j,iEl,iVar)+ & + (interp%bMatrix(i,2)*(+1.0_prec)+ & ! east (+n) + interp%bMatrix(i,1)*(-1.0_prec))/ & ! west (-n) + interp%qWeights(i)+ & + (interp%bMatrix(j,2)*(+1.0_prec)+ & ! north (+n) + interp%bMatrix(j,1)*(-1.0_prec))/ & ! south (-n) + interp%qWeights(j) + + enddo + + df%interior = abs(df%interior-0.0_prec) + + if(maxval(df%interior) <= tolerance) then + r = 0 + else + print*,"absmax error :",maxval(df%interior) + r = 1 + endif + + call f%free() + call df%free() + call interp%free() + + endfunction twopointvectordivergence_2d_constant + +endprogram test diff --git a/test/twopointvectordivergence_2d_linear.f90 b/test/twopointvectordivergence_2d_linear.f90 new file mode 100644 index 000000000..e46226946 --- /dev/null +++ b/test/twopointvectordivergence_2d_linear.f90 @@ -0,0 +1,141 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = twopointvectordivergence_2d_linear() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function twopointvectordivergence_2d_linear() result(r) + !! Verifies that the split-form divergence is exact for the field + !! V = xi1 * e_xi1 + xi2 * e_xi2 + !! on the reference element, for which the exact divergence is 2. + !! + !! The two-point arithmetic-mean flux for each computational direction is: + !! F^1(nn,i,j) = (xi1_i + xi1_nn) / 2 (xi^1 sum pairs) + !! F^2(nn,i,j) = (xi2_j + xi2_nn) / 2 (xi^2 sum pairs) + !! + !! Because the split-form sum with an arithmetic-mean two-point flux + !! reproduces the standard spectral-element derivative exactly, this test + !! verifies that the Divergence kernel correctly applies the derivative + !! matrix with the factor-of-two scaling. + use SELF_Constants + use SELF_Lagrange + use SELF_Scalar_2D + use SELF_TwoPointVector_2D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 3 + integer,parameter :: nelem = 100 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(TwoPointVector2D) :: f + type(Scalar2D) :: df + type(Lagrange),target :: interp + integer :: i,j,nn,iel,ivar + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call f%Init(interp,nvar,nelem) + call df%Init(interp,nvar,nelem) + + ! Set up two-point arithmetic-mean fluxes for V = (xi1, xi2). + ! interior(nn,i,j,iel,ivar,1) = xi^1 two-point flux for pair (i,j)-(nn,j) + ! interior(nn,i,j,iel,ivar,2) = xi^2 two-point flux for pair (i,j)-(i,nn) + do ivar = 1,nvar + do iel = 1,nelem + do j = 1,controlDegree+1 + do i = 1,controlDegree+1 + do nn = 1,controlDegree+1 + f%interior(nn,i,j,iel,ivar,1) = & + 0.5_prec*(interp%controlPoints(i)+interp%controlPoints(nn)) + f%interior(nn,i,j,iel,ivar,2) = & + 0.5_prec*(interp%controlPoints(j)+interp%controlPoints(nn)) + enddo + enddo + enddo + enddo + enddo + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%Divergence(df%interior_gpu) +#else + call f%Divergence(df%interior) +#endif + call df%UpdateHost() + + ! Add surface term. For V = (xi1, xi2), the physical flux at each face is: + ! east: f_xi1 at xi^1=+1 = cp(N+1), fbn = +cp(N+1) + ! west: f_xi1 at xi^1=-1 = cp(1), fbn = -cp(1) + ! north: f_xi2 at xi^2=+1 = cp(N+1), fbn = +cp(N+1) + ! south: f_xi2 at xi^2=-1 = cp(1), fbn = -cp(1) + do concurrent(i=1:controlDegree+1,j=1:controlDegree+1, & + iEl=1:nelem,iVar=1:nvar) + + df%interior(i,j,iEl,iVar) = df%interior(i,j,iEl,iVar)+ & + (interp%bMatrix(i,2)*(+interp%controlPoints(controlDegree+1))+ & + interp%bMatrix(i,1)*(-interp%controlPoints(1)))/ & + interp%qWeights(i)+ & + (interp%bMatrix(j,2)*(+interp%controlPoints(controlDegree+1))+ & + interp%bMatrix(j,1)*(-interp%controlPoints(1)))/ & + interp%qWeights(j) + + enddo + + print*,"absmax (nabla.F) :",maxval(df%interior) + df%interior = abs(df%interior-2.0_prec) + print*,"absmax error :",maxval(df%interior) + + if(maxval(df%interior) <= tolerance) then + r = 0 + else + r = 1 + endif + + call f%free() + call df%free() + call interp%free() + + endfunction twopointvectordivergence_2d_linear + +endprogram test diff --git a/test/twopointvectordivergence_2d_rotational.f90 b/test/twopointvectordivergence_2d_rotational.f90 new file mode 100644 index 000000000..211d99510 --- /dev/null +++ b/test/twopointvectordivergence_2d_rotational.f90 @@ -0,0 +1,137 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = twopointvectordivergence_2d_rotational() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function twopointvectordivergence_2d_rotational() result(r) + !! Verifies that the split-form divergence is identically zero for the + !! purely rotational field V = -xi2 * e_xi1 + xi1 * e_xi2. + !! + !! The exact divergence is d(-xi2)/d(xi1) + d(xi1)/d(xi2) = 0. + !! + !! Two-point arithmetic-mean fluxes (linear components, no projection error): + !! F^1(nn,i,j) = (-xi2_j + (-xi2_j))/2 = -xi2_j (independent of nn) + !! F^2(nn,i,j) = (xi1_i + xi1_i)/2 = xi1_i (independent of nn) + !! + !! Since neither flux depends on the two-point index nn, the split-form + !! sum reduces to the standard D-matrix column sum, which is zero. + use SELF_Constants + use SELF_Lagrange + use SELF_Scalar_2D + use SELF_TwoPointVector_2D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 3 + integer,parameter :: nelem = 100 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(TwoPointVector2D) :: f + type(Scalar2D) :: df + type(Lagrange),target :: interp + integer :: i,j,nn,iel,ivar + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call f%Init(interp,nvar,nelem) + call df%Init(interp,nvar,nelem) + + ! V = (-xi2, xi1): two-point fluxes are independent of nn. + ! F^1(nn,i,j) = -xi2_j (xi1-component doesn't depend on xi1) + ! F^2(nn,i,j) = xi1_i (xi2-component doesn't depend on xi2) + do ivar = 1,nvar + do iel = 1,nelem + do j = 1,controlDegree+1 + do i = 1,controlDegree+1 + do nn = 1,controlDegree+1 + f%interior(nn,i,j,iel,ivar,1) = -interp%controlPoints(j) + f%interior(nn,i,j,iel,ivar,2) = interp%controlPoints(i) + enddo + enddo + enddo + enddo + enddo + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%Divergence(df%interior_gpu) +#else + call f%Divergence(df%interior) +#endif + call df%UpdateHost() + + ! Add surface term. For V = (-xi2, xi1): + ! east: F^xi1 = -cp(j), fbn = +(-cp(j)) = -cp(j) + ! west: F^xi1 = -cp(j), fbn = -(-cp(j)) = +cp(j) + ! north: F^xi2 = +cp(i), fbn = +cp(i) + ! south: F^xi2 = +cp(i), fbn = -cp(i) + do concurrent(i=1:controlDegree+1,j=1:controlDegree+1, & + iEl=1:nelem,iVar=1:nvar) + + df%interior(i,j,iEl,iVar) = df%interior(i,j,iEl,iVar)+ & + (interp%bMatrix(i,2)*(-interp%controlPoints(j))+ & + interp%bMatrix(i,1)*(+interp%controlPoints(j)))/ & + interp%qWeights(i)+ & + (interp%bMatrix(j,2)*(+interp%controlPoints(i))+ & + interp%bMatrix(j,1)*(-interp%controlPoints(i)))/ & + interp%qWeights(j) + + enddo + + df%interior = abs(df%interior-0.0_prec) + + if(maxval(df%interior) <= tolerance) then + r = 0 + else + print*,"absmax error :",maxval(df%interior) + r = 1 + endif + + call f%free() + call df%free() + call interp%free() + + endfunction twopointvectordivergence_2d_rotational + +endprogram test diff --git a/test/twopointvectordivergence_3d_constant.f90 b/test/twopointvectordivergence_3d_constant.f90 new file mode 100644 index 000000000..6d2236273 --- /dev/null +++ b/test/twopointvectordivergence_3d_constant.f90 @@ -0,0 +1,112 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = twopointvectordivergence_3d_constant() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function twopointvectordivergence_3d_constant() result(r) + !! Verifies that the 3-D split-form divergence of a constant two-point + !! vector field is identically zero. + use SELF_Constants + use SELF_Lagrange + use SELF_Scalar_3D + use SELF_TwoPointVector_3D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 3 + integer,parameter :: nelem = 64 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(TwoPointVector3D) :: f + type(Scalar3D) :: df + type(Lagrange),target :: interp + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call f%Init(interp,nvar,nelem) + call df%Init(interp,nvar,nelem) + + f%interior = 1.0_prec + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%Divergence(df%interior_gpu) +#else + call f%Divergence(df%interior) +#endif + call df%UpdateHost() + + ! Add surface term for each of three coordinate directions + block + integer :: i,j,k,iEl,iVar + do concurrent(i=1:controlDegree+1,j=1:controlDegree+1, & + k=1:controlDegree+1,iEl=1:nelem,iVar=1:nvar) + + df%interior(i,j,k,iEl,iVar) = df%interior(i,j,k,iEl,iVar)+ & + (interp%bMatrix(i,2)*(+1.0_prec)+ & + interp%bMatrix(i,1)*(-1.0_prec))/interp%qWeights(i)+ & + (interp%bMatrix(j,2)*(+1.0_prec)+ & + interp%bMatrix(j,1)*(-1.0_prec))/interp%qWeights(j)+ & + (interp%bMatrix(k,2)*(+1.0_prec)+ & + interp%bMatrix(k,1)*(-1.0_prec))/interp%qWeights(k) + + enddo + endblock + + df%interior = abs(df%interior-0.0_prec) + + if(maxval(df%interior) <= tolerance) then + r = 0 + else + print*,"absmax error :",maxval(df%interior) + r = 1 + endif + + call f%free() + call df%free() + call interp%free() + + endfunction twopointvectordivergence_3d_constant + +endprogram test diff --git a/test/twopointvectordivergence_3d_linear.f90 b/test/twopointvectordivergence_3d_linear.f90 new file mode 100644 index 000000000..8843f0c29 --- /dev/null +++ b/test/twopointvectordivergence_3d_linear.f90 @@ -0,0 +1,135 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program test + + implicit none + integer :: exit_code + + exit_code = twopointvectordivergence_3d_linear() + if(exit_code /= 0) then + stop exit_code + endif + +contains + + integer function twopointvectordivergence_3d_linear() result(r) + !! Verifies that the 3-D split-form divergence is exact for the field + !! V = xi1 * e_xi1 + xi2 * e_xi2 + xi3 * e_xi3 + !! on the reference element, for which the exact divergence is 3. + use SELF_Constants + use SELF_Lagrange + use SELF_Scalar_3D + use SELF_TwoPointVector_3D + + implicit none + + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + integer,parameter :: nvar = 1 + integer,parameter :: nelem = 64 +#ifdef DOUBLE_PRECISION + real(prec),parameter :: tolerance = 10.0_prec**(-7) +#else + real(prec),parameter :: tolerance = 10.0_prec**(-3) +#endif + type(TwoPointVector3D) :: f + type(Scalar3D) :: df + type(Lagrange),target :: interp + integer :: i,j,k,nn,iel,ivar + + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS_LOBATTO, & + M=targetDegree, & + targetNodeType=UNIFORM) + + call f%Init(interp,nvar,nelem) + call df%Init(interp,nvar,nelem) + + ! Two-point arithmetic-mean fluxes for V = (xi1, xi2, xi3). + ! F^r(nn,i,j,k) = (coord_r_at_node + coord_r_at_nn) / 2 + do ivar = 1,nvar + do iel = 1,nelem + do k = 1,controlDegree+1 + do j = 1,controlDegree+1 + do i = 1,controlDegree+1 + do nn = 1,controlDegree+1 + f%interior(nn,i,j,k,iel,ivar,1) = & + 0.5_prec*(interp%controlPoints(i)+interp%controlPoints(nn)) + f%interior(nn,i,j,k,iel,ivar,2) = & + 0.5_prec*(interp%controlPoints(j)+interp%controlPoints(nn)) + f%interior(nn,i,j,k,iel,ivar,3) = & + 0.5_prec*(interp%controlPoints(k)+interp%controlPoints(nn)) + enddo + enddo + enddo + enddo + enddo + enddo + call f%UpdateDevice() + +#ifdef ENABLE_GPU + call f%Divergence(df%interior_gpu) +#else + call f%Divergence(df%interior) +#endif + call df%UpdateHost() + + ! Add surface term for V = (xi1, xi2, xi3) + block + integer :: i,j,k,iEl,iVar + integer :: Np1 + Np1 = controlDegree+1 + do concurrent(i=1:Np1,j=1:Np1,k=1:Np1,iEl=1:nelem,iVar=1:nvar) + + df%interior(i,j,k,iEl,iVar) = df%interior(i,j,k,iEl,iVar)+ & + (interp%bMatrix(i,2)*(+interp%controlPoints(Np1))+ & + interp%bMatrix(i,1)*(-interp%controlPoints(1)))/interp%qWeights(i)+ & + (interp%bMatrix(j,2)*(+interp%controlPoints(Np1))+ & + interp%bMatrix(j,1)*(-interp%controlPoints(1)))/interp%qWeights(j)+ & + (interp%bMatrix(k,2)*(+interp%controlPoints(Np1))+ & + interp%bMatrix(k,1)*(-interp%controlPoints(1)))/interp%qWeights(k) + + enddo + endblock + + print*,"absmax (nabla.F) :",maxval(df%interior) + df%interior = abs(df%interior-3.0_prec) + print*,"absmax error :",maxval(df%interior) + + if(maxval(df%interior) <= tolerance) then + r = 0 + else + r = 1 + endif + + call f%free() + call df%free() + call interp%free() + + endfunction twopointvectordivergence_3d_linear + +endprogram test