! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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 use SELF_BoundaryConditions 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 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. !! 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 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) !! 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(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) = & m%solution%boundary(i,j,iEl,1:m%nvar) enddo enddo endselect endsubroutine hbc2d_NoNormalFlow_ECAdvection2D endmodule SELF_ECAdvection2D_t