! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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 use SELF_mesh use SELF_BoundaryConditions implicit none type,extends(ECDGModel2D) :: 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 procedure :: AdditionalInit => AdditionalInit_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. !! Uses the per-face spectral radius lambda = |u.n|, matching upwind !! (Godunov) for linear advection: tangential faces (u.n=0) contribute !! zero flux, avoiding spurious dissipation across element interfaces !! whose face normal is perpendicular to the velocity. 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 = abs(un) 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 subroutine AdditionalInit_ECAdvection2D_t(this) implicit none class(ECAdvection2D_t),intent(inout) :: this ! Local procedure(SELF_bcMethod),pointer :: bcfunc bcfunc => hbc2d_NoNormalFlow_ECAdvection2D call this%hyperbolicBCs%RegisterBoundaryCondition( & SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) endsubroutine AdditionalInit_ECAdvection2D_t subroutine hbc2d_NoNormalFlow_ECAdvection2D(bc,mymodel) !! Prescribed-zero boundary state for the tracer. !! !! "No-normal-flow" only has physical meaning when the velocity field !! is a prognostic variable that can be reflected at the wall. For a !! passive tracer with externally-prescribed velocity the appropriate !! external state is identically zero everywhere on the boundary, !! which combined with the LLF (= upwind for linear advection) Riemann !! flux gives correct inflow (zero injection) and outflow (upwind sL) !! behaviour. class(BoundaryCondition),intent(in) :: bc class(Model),intent(inout) :: mymodel ! Local integer :: n,i,iEl,j select type(m => mymodel) class is(ECAdvection2D_t) do n = 1,bc%nBoundaries iEl = bc%elements(n) j = bc%sides(n) do i = 1,m%solution%interp%N+1 m%solution%extBoundary(i,j,iEl,1:m%nvar) = 0.0_prec enddo enddo endselect endsubroutine hbc2d_NoNormalFlow_ECAdvection2D endmodule SELF_ECAdvection2D_t