! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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_DGModel2D_t use SELF_SupportRoutines use SELF_Metadata use SELF_Geometry_2D use SELF_Mesh_2D use SELF_MappedScalar_2D use SELF_MappedVector_2D use SELF_HDF5 use HDF5 use FEQParse use SELF_Model use SELF_BoundaryConditions use SELF_TransferPlan_2D implicit none type,extends(Model) :: DGModel2D_t type(MappedScalar2D) :: solution type(MappedVector2D) :: solutionGradient type(MappedVector2D) :: flux type(MappedScalar2D) :: source type(MappedScalar2D) :: fluxDivergence type(MappedScalar2D) :: dSdt type(MappedScalar2D) :: workSol type(Mesh2D),pointer :: mesh => null() type(SEMQuad),pointer :: geometry => null() type(BoundaryConditionList) :: hyperbolicBCs type(BoundaryConditionList) :: parabolicBCs !! Pre-regrid copy of the solution, held between StageSolutionForTransfer and !! ApplyTransferPlan so that Regrid is free to release the storage it was read from. The !! base implementation stages on the host; the GPU backend overrides both procedures and !! stages device-side instead, leaving this unallocated. real(prec),allocatable :: transferStage(:,:,:,:) contains procedure :: Init => Init_DGModel2D_t procedure :: SetMetadata => SetMetadata_DGModel2D_t procedure :: Free => Free_DGModel2D_t procedure :: Regrid => Regrid_DGModel2D_t procedure :: MapBoundaryConditions => MapBoundaryConditions_DGModel2D_t procedure :: StageSolutionForTransfer => StageSolutionForTransfer_DGModel2D_t procedure :: ApplyTransferPlan => ApplyTransferPlan_DGModel2D_t procedure :: CalculateEntropy => CalculateEntropy_DGModel2D_t procedure :: BoundaryFlux => BoundaryFlux_DGModel2D_t procedure :: FluxMethod => fluxmethod_DGModel2D_t procedure :: SourceMethod => sourcemethod_DGModel2D_t procedure :: SetBoundaryCondition => setboundarycondition_DGModel2D_t procedure :: SetGradientBoundaryCondition => setgradientboundarycondition_DGModel2D_t procedure :: ReportMetrics => ReportMetrics_DGModel2D_t procedure :: UpdateSolution => UpdateSolution_DGModel2D_t procedure :: UpdateGRK2 => UpdateGRK2_DGModel2D_t procedure :: UpdateGRK3 => UpdateGRK3_DGModel2D_t procedure :: UpdateGRK4 => UpdateGRK4_DGModel2D_t procedure :: CalculateSolutionGradient => CalculateSolutionGradient_DGModel2D_t procedure :: CalculateTendency => CalculateTendency_DGModel2D_t generic :: SetSolution => SetSolutionFromChar_DGModel2D_t, & SetSolutionFromEqn_DGModel2D_t procedure,private :: SetSolutionFromChar_DGModel2D_t procedure,private :: SetSolutionFromEqn_DGModel2D_t procedure :: ReadModel => Read_DGModel2D_t procedure :: WriteModel => Write_DGModel2D_t procedure :: WriteTecplot => WriteTecplot_DGModel2D_t endtype DGModel2D_t contains subroutine Init_DGModel2D_t(this,mesh,geometry) implicit none class(DGModel2D_t),intent(out) :: this type(Mesh2D),intent(in),target :: mesh type(SEMQuad),intent(in),target :: geometry ! Local this%mesh => mesh this%geometry => geometry call this%SetNumberOfVariables() ! Default the number of time-stepped variables to nvar. Models that carry ! auxiliary/diagnostic variables may set this%nstepped < nvar inside ! SetNumberOfVariables to exclude the trailing variables from time integration. if(this%nstepped <= 0 .or. this%nstepped > this%nvar) this%nstepped = this%nvar call this%solution%Init(geometry%x%interp,this%nvar,this%mesh%nElem) call this%workSol%Init(geometry%x%interp,this%nvar,this%mesh%nElem) call this%dSdt%Init(geometry%x%interp,this%nvar,this%mesh%nElem) call this%solutionGradient%Init(geometry%x%interp,this%nvar,this%mesh%nElem) call this%flux%Init(geometry%x%interp,this%nvar,this%mesh%nElem) call this%source%Init(geometry%x%interp,this%nvar,this%mesh%nElem) call this%fluxDivergence%Init(geometry%x%interp,this%nvar,this%mesh%nElem) call this%solution%AssociateGeometry(geometry) call this%solutionGradient%AssociateGeometry(geometry) call this%flux%AssociateGeometry(geometry) call this%fluxDivergence%AssociateGeometry(geometry) call this%hyperbolicBCs%Init() call this%parabolicBCs%Init() call this%AdditionalInit() call this%MapBoundaryConditions() call this%SetMetadata() endsubroutine Init_DGModel2D_t subroutine SetMetadata_DGModel2D_t(this) implicit none class(DGModel2D_t),intent(inout) :: this ! Local integer :: ivar character(LEN=3) :: ivarChar character(LEN=25) :: varname do ivar = 1,this%nvar write(ivarChar,'(I3.3)') ivar varname = "solution"//trim(ivarChar) call this%solution%SetName(ivar,varname) call this%solution%SetUnits(ivar,"[null]") enddo endsubroutine SetMetadata_DGModel2D_t subroutine Free_DGModel2D_t(this) implicit none class(DGModel2D_t),intent(inout) :: this call this%solution%Free() call this%workSol%Free() call this%dSdt%Free() call this%solutionGradient%Free() call this%flux%Free() call this%source%Free() call this%fluxDivergence%Free() call this%hyperbolicBCs%Free() call this%parabolicBCs%Free() call this%AdditionalFree() if(allocated(this%transferStage)) deallocate(this%transferStage) endsubroutine Free_DGModel2D_t subroutine StageSolutionForTransfer_DGModel2D_t(this) !! Preserve the current solution ahead of a regrid, so that Regrid may release the storage !! it lives in. Pair with ApplyTransferPlan, which consumes the staged copy: !! !! call model%StageSolutionForTransfer() !! call model%Regrid(newMesh,newGeom) !! call model%ApplyTransferPlan(plan,interp,eFirst,eLast) !! !! This base implementation stages on the host, which on a GPU build means a !! device-to-host copy of the whole field; the GPU backend overrides it with a !! device-to-device copy and no host traffic (Stage 6a). implicit none class(DGModel2D_t),intent(inout) :: this ! Local integer :: Np,nEl Np = this%solution%interp%N+1 nEl = this%solution%nElem call this%solution%UpdateHost() if(allocated(this%transferStage)) deallocate(this%transferStage) allocate(this%transferStage(1:Np,1:Np,1:nEl,1:this%nvar)) this%transferStage(1:Np,1:Np,1:nEl,1:this%nvar) = & this%solution%interior(1:Np,1:Np,1:nEl,1:this%nvar) endsubroutine StageSolutionForTransfer_DGModel2D_t subroutine ApplyTransferPlan_DGModel2D_t(this,plan,interp,eFirst,eLast,uGlobal) !! Transfer the staged pre-regrid solution onto the regridded mesh through plan, filling the !! rank-local element range [eFirst,eLast] of the new solution. !! !! uGlobal is optional and supplies the GLOBAL old field when the caller has already !! assembled one (the multi-rank allgather path); when absent the locally staged copy from !! StageSolutionForTransfer is used, which is the whole field on a single rank. !! !! This base implementation runs the portable host transfer and uploads the result; the GPU !! backend overrides it to run the transfer on the device with no host traffic. implicit none class(DGModel2D_t),intent(inout) :: this ! target: the GPU override takes c_loc of the plan's arrays to upload them, which requires ! the POINTER or TARGET attribute. Declared here too so the override's characteristics match. type(TransferPlan2D),intent(in),target :: plan type(Lagrange),intent(in) :: interp integer,intent(in) :: eFirst integer,intent(in) :: eLast real(prec),intent(in),optional :: uGlobal(:,:,:,:) if(present(uGlobal)) then call ApplyTransferPlanRange(plan,interp,this%nvar,uGlobal,eFirst,eLast, & this%solution%interior) else if(.not. allocated(this%transferStage)) then print*,__FILE__,':',__LINE__, & ' : Error : ApplyTransferPlan called without a staged solution.' stop 1 endif call ApplyTransferPlanRange(plan,interp,this%nvar,this%transferStage,eFirst,eLast, & this%solution%interior) endif call this%solution%UpdateDevice() if(allocated(this%transferStage)) deallocate(this%transferStage) endsubroutine ApplyTransferPlan_DGModel2D_t subroutine Regrid_DGModel2D_t(this,mesh,geometry) !! Rebind a live model to a new mesh/geometry pair (AMR regrid). The mesh-sized solution !! storage is reallocated and the boundary-condition registrations and maps are rebuilt !! for the new mesh, while everything that is not mesh-sized is preserved: the time state !! (t, dt, entropy, IO counter), the time-integrator selection, configuration flags, and !! any model-specific parameters (Init is intent(out) and would reset all of these). !! nvar/nstepped are unchanged - the model solves the same equations on a new mesh. !! !! The solution interior is left UNINITIALIZED: the caller transfers the solution from the !! previous mesh (e.g. ApplyTransferPlan on a BuildTransferPlan mapping) and then calls !! solution%UpdateDevice. Regrid runs once per adaptation epoch, between time steps; it is !! not a per-step hot path. implicit none class(DGModel2D_t),intent(inout) :: this type(Mesh2D),intent(in),target :: mesh type(SEMQuad),intent(in),target :: geometry if(.not. associated(this%mesh)) then print*,__FILE__,':',__LINE__, & ' : Error : Regrid called on a model that has not been initialized.' stop 1 endif ! Free everything sized by the old mesh, mirroring Free (AdditionalFree releases any ! model-specific mesh-sized state so AdditionalInit can rebuild it below). ! Boundary-condition registrations are rebuilt because the boundary side set changes with the ! mesh. The mesh-sized fields are NOT freed: they are resized in place below (AMR Stage 6b), ! which reuses their host pools and device buffers whenever the new element count fits. call this%hyperbolicBCs%Free() call this%parabolicBCs%Free() call this%AdditionalFree() ! Rebuild on the new mesh, mirroring the mesh-sized portion of Init. this%mesh => mesh this%geometry => geometry ! Resize rather than Free + Init. Init is intent(out), so it would reset the whole object, ! reallocate every array, zero it, reconstruct the equation parsers and - on GPU builds - ! upload the zeros, all of which the adaptive loop then discards. Profiling attributed over ! half of an adaptation to exactly that cycle. call this%solution%Resize(geometry%x%interp,this%nvar,this%mesh%nElem) call this%workSol%Resize(geometry%x%interp,this%nvar,this%mesh%nElem) call this%dSdt%Resize(geometry%x%interp,this%nvar,this%mesh%nElem) call this%solutionGradient%Resize(geometry%x%interp,this%nvar,this%mesh%nElem) call this%flux%Resize(geometry%x%interp,this%nvar,this%mesh%nElem) call this%source%Resize(geometry%x%interp,this%nvar,this%mesh%nElem) call this%fluxDivergence%Resize(geometry%x%interp,this%nvar,this%mesh%nElem) call this%solution%AssociateGeometry(geometry) call this%solutionGradient%AssociateGeometry(geometry) call this%flux%AssociateGeometry(geometry) call this%fluxDivergence%AssociateGeometry(geometry) call this%hyperbolicBCs%Init() call this%parabolicBCs%Init() call this%AdditionalInit() call this%MapBoundaryConditions() call this%SetMetadata() endsubroutine Regrid_DGModel2D_t subroutine ReportMetrics_DGModel2D_t(this) !! Base method for reporting the entropy of a model !! to stdout. Only override this procedure if additional !! reporting is needed. Alternatively, if you think !! additional reporting would be valuable for all models, !! open a pull request with modifications to this base !! method. implicit none class(DGModel2D_t),intent(inout) :: this ! Local character(len=20) :: modelTime character(len=20) :: minv,maxv character(len=:),allocatable :: str integer :: ivar ! Copy the time and entropy to a string write(modelTime,"(ES16.7E3)") this%t do ivar = 1,this%nvar write(maxv,"(ES16.7E3)") maxval(this%solution%interior(:,:,:,ivar)) write(minv,"(ES16.7E3)") minval(this%solution%interior(:,:,:,ivar)) ! Write the output to STDOUT open(output_unit,ENCODING='utf-8') write(output_unit,'(1x, A," : ")',ADVANCE='no') __FILE__ str = 'tᵢ ='//trim(modelTime) write(output_unit,'(A)',ADVANCE='no') str str = ' | min('//trim(this%solution%meta(ivar)%name)// & '), max('//trim(this%solution%meta(ivar)%name)//') = '// & minv//" , "//maxv write(output_unit,'(A)',ADVANCE='yes') str enddo call this%ReportUserMetrics() endsubroutine ReportMetrics_DGModel2D_t subroutine SetSolutionFromEqn_DGModel2D_t(this,eqn) implicit none class(DGModel2D_t),intent(inout) :: this type(EquationParser),intent(in) :: eqn(1:this%solution%nVar) ! Local integer :: iVar ! Copy the equation parser do iVar = 1,this%solution%nVar call this%solution%SetEquation(ivar,eqn(iVar)%equation) enddo call this%solution%SetInteriorFromEquation(this%geometry,this%t) call this%solution%BoundaryInterp() endsubroutine SetSolutionFromEqn_DGModel2D_t subroutine SetSolutionFromChar_DGModel2D_t(this,eqnChar) implicit none class(DGModel2D_t),intent(inout) :: this character(*),intent(in) :: eqnChar(1:this%solution%nVar) ! Local integer :: iVar do iVar = 1,this%solution%nVar call this%solution%SetEquation(ivar,trim(eqnChar(iVar))) enddo call this%solution%SetInteriorFromEquation(this%geometry,this%t) call this%solution%BoundaryInterp() endsubroutine SetSolutionFromChar_DGModel2D_t subroutine UpdateSolution_DGModel2D_t(this,dt) !! Computes a solution update as , where dt is either provided through the interface !! or taken as the Model's stored time step size (model % dt) implicit none class(DGModel2D_t),intent(inout) :: this real(prec),optional,intent(in) :: dt ! Local real(prec) :: dtLoc integer :: i,j,iEl,iVar if(present(dt)) then dtLoc = dt else dtLoc = this%dt endif do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & iel=1:this%mesh%nElem,ivar=1:this%nstepped) this%solution%interior(i,j,iEl,iVar) = & this%solution%interior(i,j,iEl,iVar)+ & dtLoc*this%dSdt%interior(i,j,iEl,iVar) enddo endsubroutine UpdateSolution_DGModel2D_t subroutine UpdateGRK2_DGModel2D_t(this,m) implicit none class(DGModel2D_t),intent(inout) :: this integer,intent(in) :: m ! Local integer :: i,j,iEl,iVar do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & iel=1:this%mesh%nElem,ivar=1:this%nstepped) this%workSol%interior(i,j,iEl,iVar) = rk2_a(m)* & this%workSol%interior(i,j,iEl,iVar)+ & this%dSdt%interior(i,j,iEl,iVar) this%solution%interior(i,j,iEl,iVar) = & this%solution%interior(i,j,iEl,iVar)+ & rk2_g(m)*this%dt*this%workSol%interior(i,j,iEl,iVar) enddo endsubroutine UpdateGRK2_DGModel2D_t subroutine UpdateGRK3_DGModel2D_t(this,m) implicit none class(DGModel2D_t),intent(inout) :: this integer,intent(in) :: m ! Local integer :: i,j,iEl,iVar do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & iel=1:this%mesh%nElem,ivar=1:this%nstepped) this%workSol%interior(i,j,iEl,iVar) = rk3_a(m)* & this%workSol%interior(i,j,iEl,iVar)+ & this%dSdt%interior(i,j,iEl,iVar) this%solution%interior(i,j,iEl,iVar) = & this%solution%interior(i,j,iEl,iVar)+ & rk3_g(m)*this%dt*this%workSol%interior(i,j,iEl,iVar) enddo endsubroutine UpdateGRK3_DGModel2D_t subroutine UpdateGRK4_DGModel2D_t(this,m) implicit none class(DGModel2D_t),intent(inout) :: this integer,intent(in) :: m ! Local integer :: i,j,iEl,iVar do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & iel=1:this%mesh%nElem,ivar=1:this%nstepped) this%workSol%interior(i,j,iEl,iVar) = rk4_a(m)* & this%workSol%interior(i,j,iEl,iVar)+ & this%dSdt%interior(i,j,iEl,iVar) this%solution%interior(i,j,iEl,iVar) = & this%solution%interior(i,j,iEl,iVar)+ & rk4_g(m)*this%dt*this%workSol%interior(i,j,iEl,iVar) enddo endsubroutine UpdateGRK4_DGModel2D_t subroutine CalculateSolutionGradient_DGModel2D_t(this) implicit none class(DGModel2D_t),intent(inout) :: this call this%solution%AverageSides() call this%solution%MappedDGGradient(this%solutionGradient%interior) ! interpolate the solutiongradient to the element boundaries call this%solutionGradient%BoundaryInterp() ! perform the side exchange to populate the ! solutionGradient % extBoundary attribute call this%solutionGradient%SideExchange(this%mesh) ! populate the solutionGradient % extBoundary attribute on ! nonconforming (mortar) interfaces if(this%mesh%nMortars > 0) then call this%solutionGradient%MortarExchange(this%mesh) endif endsubroutine CalculateSolutionGradient_DGModel2D_t subroutine CalculateEntropy_DGModel2D_t(this) implicit none class(DGModel2D_t),intent(inout) :: this ! Local integer :: iel,i,j,ierror real(prec) :: e,jac real(prec) :: s(1:this%nvar) e = 0.0_prec do iel = 1,this%geometry%nelem do j = 1,this%solution%interp%N+1 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* & this%solution%interp%qWeights(i)* & this%solution%interp%qWeights(j) enddo enddo enddo if(this%mesh%decomp%mpiEnabled) then call mpi_allreduce(e, & this%entropy, & 1, & this%mesh%decomp%mpiPrec, & MPI_SUM, & this%mesh%decomp%mpiComm, & iError) else this%entropy = e endif endsubroutine CalculateEntropy_DGModel2D_t subroutine fluxmethod_DGModel2D_t(this) implicit none class(DGModel2D_t),intent(inout) :: this ! Local integer :: iel integer :: i integer :: j real(prec) :: s(1:this%nvar),dsdx(1:this%nvar,1:2) do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & iel=1:this%mesh%nElem) s = this%solution%interior(i,j,iel,1:this%nvar) dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) this%flux%interior(i,j,iel,1:this%nvar,1:2) = this%flux2d(s,dsdx) enddo endsubroutine fluxmethod_DGModel2D_t subroutine BoundaryFlux_DGModel2D_t(this) ! this method uses an linear upwind solver for the ! advective flux and the bassi-rebay method for the ! diffusive fluxes implicit none class(DGModel2D_t),intent(inout) :: this ! Local integer :: iel integer :: j integer :: i real(prec) :: sL(1:this%nvar),sR(1:this%nvar) real(prec) :: dsdx(1:this%nvar,1:2) real(prec) :: nhat(1:2),nmag do concurrent(i=1:this%solution%N+1,j=1:4, & iel=1:this%mesh%nElem) ! Get the boundary normals on cell edges from the mesh geometry nhat = this%geometry%nHat%boundary(i,j,iEl,1,1:2) sL = this%solution%boundary(i,j,iel,1:this%nvar) ! interior solution sR = this%solution%extboundary(i,j,iel,1:this%nvar) ! exterior solution dsdx = this%solutiongradient%avgboundary(i,j,iel,1:this%nvar,1:2) nmag = this%geometry%nScale%boundary(i,j,iEl,1) this%flux%boundaryNormal(i,j,iEl,1:this%nvar) = this%riemannflux2d(sL,sR,dsdx,nhat)*nmag enddo endsubroutine BoundaryFlux_DGModel2D_t subroutine sourcemethod_DGModel2D_t(this) implicit none class(DGModel2D_t),intent(inout) :: this ! Local integer :: iel integer :: i integer :: j real(prec) :: s(1:this%nvar),dsdx(1:this%nvar,1:2) do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, & iel=1:this%mesh%nElem) s = this%solution%interior(i,j,iel,1:this%nvar) dsdx = this%solutionGradient%interior(i,j,iel,1:this%nvar,1:2) this%source%interior(i,j,iel,1:this%nvar) = this%source2d(s,dsdx) enddo endsubroutine sourcemethod_DGModel2D_t subroutine MapBoundaryConditions_DGModel2D_t(this) !! Scan the mesh sideInfo and populate the elements/sides !! arrays for each registered boundary condition. implicit none class(DGModel2D_t),intent(inout) :: this ! Local type(BoundaryCondition),pointer :: bc integer :: iEl,j,e2,bcid integer :: count,n integer,allocatable :: elems(:),sds(:) ! Map hyperbolic BCs bc => this%hyperbolicBCs%head do while(associated(bc)) ! Pass 1: count boundary faces for this bcid count = 0 do iEl = 1,this%mesh%nElem do j = 1,4 e2 = this%mesh%sideInfo(3,j,iEl) bcid = this%mesh%sideInfo(5,j,iEl) if(e2 == 0 .and. bcid == bc%bcid) count = count+1 enddo enddo if(count > 0) then ! Pass 2: fill element/side arrays allocate(elems(count),sds(count)) n = 0 do iEl = 1,this%mesh%nElem do j = 1,4 e2 = this%mesh%sideInfo(3,j,iEl) bcid = this%mesh%sideInfo(5,j,iEl) if(e2 == 0 .and. bcid == bc%bcid) then n = n+1 elems(n) = iEl sds(n) = j endif enddo enddo call this%hyperbolicBCs%PopulateBoundaries(bc%bcid,count,elems,sds) deallocate(elems,sds) endif bc => bc%next enddo ! Map parabolic BCs bc => this%parabolicBCs%head do while(associated(bc)) count = 0 do iEl = 1,this%mesh%nElem do j = 1,4 e2 = this%mesh%sideInfo(3,j,iEl) bcid = this%mesh%sideInfo(5,j,iEl) if(e2 == 0 .and. bcid == bc%bcid) count = count+1 enddo enddo if(count > 0) then allocate(elems(count),sds(count)) n = 0 do iEl = 1,this%mesh%nElem do j = 1,4 e2 = this%mesh%sideInfo(3,j,iEl) bcid = this%mesh%sideInfo(5,j,iEl) if(e2 == 0 .and. bcid == bc%bcid) then n = n+1 elems(n) = iEl sds(n) = j endif enddo enddo call this%parabolicBCs%PopulateBoundaries(bc%bcid,count,elems,sds) deallocate(elems,sds) endif bc => bc%next enddo endsubroutine MapBoundaryConditions_DGModel2D_t subroutine setboundarycondition_DGModel2D_t(this) !! Apply registered boundary conditions for the solution. !! Each boundary condition method loops over its own !! boundary faces. implicit none class(DGModel2D_t),intent(inout) :: this ! Local type(BoundaryCondition),pointer :: bc procedure(SELF_bcMethod),pointer :: apply_bc bc => this%hyperbolicBCs%head do while(associated(bc)) apply_bc => bc%bcMethod call apply_bc(bc,this) bc => bc%next enddo endsubroutine setboundarycondition_DGModel2D_t subroutine setgradientboundarycondition_DGModel2D_t(this) !! Apply registered boundary conditions for the solution gradient. !! Each boundary condition method loops over its own !! boundary faces. implicit none class(DGModel2D_t),intent(inout) :: this ! Local type(BoundaryCondition),pointer :: bc procedure(SELF_bcMethod),pointer :: apply_bc bc => this%parabolicBCs%head do while(associated(bc)) apply_bc => bc%bcMethod call apply_bc(bc,this) bc => bc%next enddo endsubroutine setgradientboundarycondition_DGModel2D_t subroutine CalculateTendency_DGModel2D_t(this) implicit none class(DGModel2D_t),intent(inout) :: this ! Local integer :: i,j,iEl,iVar call this%solution%BoundaryInterp() call this%solution%SideExchange(this%mesh) ! populate the solution % extBoundary attribute on nonconforming ! (mortar) interfaces if(this%mesh%nMortars > 0) then call this%solution%MortarExchange(this%mesh) endif call this%PreTendencyHook() ! User-supplied call this%SetBoundaryCondition() ! User-supplied if(this%gradient_enabled) then call this%CalculateSolutionGradient() call this%SetGradientBoundaryCondition() ! User-supplied call this%solutionGradient%AverageSides() endif call this%SourceMethod() ! User supplied call this%BoundaryFlux() ! User supplied ! On mortar interfaces, replace the big side's surface-flux integrand with the ! projection of the small sides' integrands so that the interface is conservative if(this%mesh%nMortars > 0) then call this%flux%MortarFluxCollect(this%mesh) endif call this%FluxMethod() ! User supplied call this%flux%MappedDGDivergence(this%fluxDivergence%interior) 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_DGModel2D_t subroutine Write_DGModel2D_t(this,fileName) implicit none class(DGModel2D_t),intent(inout) :: this character(*),optional,intent(in) :: fileName ! Local integer(HID_T) :: fileId character(LEN=self_FileNameLength) :: pickupFile character(13) :: timeStampString if(present(filename)) then pickupFile = filename else write(timeStampString,'(I13.13)') this%ioIterate pickupFile = 'solution.'//timeStampString//'.h5' endif print*,__FILE__//" : Writing pickup file : "//trim(pickupFile) call this%solution%UpdateHost() if(this%mesh%decomp%mpiEnabled) then call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId,this%mesh%decomp%mpiComm) ! Write the interpolant to the file call this%solution%interp%WriteHDF5(fileId) ! In this section, we write the solution and geometry on the control (quadrature) grid ! which can be used for model pickup runs or post-processing ! Write the model state to file call CreateGroup_HDF5(fileId,'/controlgrid') print*," offset, nglobal_elem : ",this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1),this%mesh%decomp%nElem call this%solution%WriteHDF5(fileId,'/controlgrid/solution', & this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1),this%mesh%decomp%nElem) ! Write the geometry to file call this%geometry%x%WriteHDF5(fileId,'/controlgrid/geometry', & this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1),this%mesh%decomp%nElem) ! -- END : writing solution on control grid -- ! call Close_HDF5(fileId) else call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId) ! Write the interpolant to the file call this%solution%interp%WriteHDF5(fileId) ! In this section, we write the solution and geometry on the control (quadrature) grid ! which can be used for model pickup runs or post-processing ! Write the model state to file call CreateGroup_HDF5(fileId,'/controlgrid') call this%solution%WriteHDF5(fileId,'/controlgrid/solution') ! Write the geometry to file call this%geometry%x%WriteHDF5(fileId,'/controlgrid/geometry') ! -- END : writing solution on control grid -- ! call Close_HDF5(fileId) endif endsubroutine Write_DGModel2D_t subroutine Read_DGModel2D_t(this,fileName) implicit none class(DGModel2D_t),intent(inout) :: this character(*),intent(in) :: fileName ! Local integer(HID_T) :: fileId integer(HID_T) :: solOffset(1:3) integer :: firstElem integer :: ivar if(this%mesh%decomp%mpiEnabled) then call Open_HDF5(fileName,H5F_ACC_RDWR_F,fileId, & this%mesh%decomp%mpiComm) else call Open_HDF5(fileName,H5F_ACC_RDWR_F,fileId) endif if(this%mesh%decomp%mpiEnabled) then firstElem = this%mesh%decomp%offsetElem(this%mesh%decomp%rankId+1) solOffset(1:3) = (/0,0,firstElem/) do ivar = 1,this%solution%nvar call ReadArray_HDF5(fileId, & '/controlgrid/solution/'//trim(this%solution%meta(ivar)%name), & this%solution%interior(:,:,:,ivar),solOffset) enddo else do ivar = 1,this%solution%nvar call ReadArray_HDF5(fileId, & '/controlgrid/solution/'//trim(this%solution%meta(ivar)%name), & this%solution%interior(:,:,:,ivar)) enddo endif call Close_HDF5(fileId) endsubroutine Read_DGModel2D_t subroutine WriteTecplot_DGModel2D_t(this,filename) implicit none class(DGModel2D_t),intent(inout) :: this character(*),intent(in),optional :: filename ! Local character(8) :: zoneID integer :: fUnit integer :: iEl,i,j,iVar character(LEN=self_FileNameLength) :: tecFile character(LEN=self_TecplotHeaderLength) :: tecHeader character(LEN=self_FormatLength) :: fmat character(13) :: timeStampString character(5) :: rankString type(Scalar2D) :: solution type(Scalar2D) :: dsdt type(Vector2D) :: solutionGradient type(Vector2D) :: x type(Lagrange),target :: interp if(present(filename)) then tecFile = filename else write(timeStampString,'(I13.13)') this%ioIterate if(this%mesh%decomp%mpiEnabled) then write(rankString,'(I5.5)') this%mesh%decomp%rankId tecFile = 'solution.'//rankString//'.'//timeStampString//'.tec' else tecFile = 'solution.'//timeStampString//'.tec' endif endif ! Create an interpolant for the uniform grid call interp%Init(this%solution%interp%M, & this%solution%interp%targetNodeType, & this%solution%interp%N, & this%solution%interp%controlNodeType) call solution%Init(interp, & this%solution%nVar,this%solution%nElem) call dsdt%Init(interp, & this%solution%nVar,this%solution%nElem) call solutionGradient%Init(interp, & this%solution%nVar,this%solution%nElem) call x%Init(interp,1,this%solution%nElem) call this%solution%UpdateHost() call this%solutionGradient%UpdateHost() call this%dsdt%UpdateHost() ! Map the mesh positions to the target grid call this%geometry%x%GridInterp(x%interior) ! Map the solution to the target grid call this%solution%GridInterp(solution%interior) call this%dsdt%GridInterp(dsdt%interior) ! Map the solution to the target grid call this%solutionGradient%GridInterp(solutionGradient%interior) open(UNIT=NEWUNIT(fUnit), & FILE=trim(tecFile), & FORM='formatted', & STATUS='replace') tecHeader = 'VARIABLES = "X", "Y"' do iVar = 1,this%solution%nVar tecHeader = trim(tecHeader)//', "'//trim(this%solution%meta(iVar)%name)//'"' enddo do iVar = 1,this%solution%nVar tecHeader = trim(tecHeader)//', "d/dx('//trim(this%solution%meta(iVar)%name)//')"' enddo do iVar = 1,this%solution%nVar tecHeader = trim(tecHeader)//', "d/dy('//trim(this%solution%meta(iVar)%name)//')"' enddo do iVar = 1,this%solution%nVar tecHeader = trim(tecHeader)//', "d/dt('//trim(this%solution%meta(iVar)%name)//')"' enddo write(fUnit,*) trim(tecHeader) ! Create format statement write(fmat,*) 4*this%solution%nvar+2 fmat = '('//trim(fmat)//'(ES16.7E3,1x))' do iEl = 1,this%solution%nElem ! TO DO :: Get the global element ID write(zoneID,'(I8.8)') iEl write(fUnit,*) 'ZONE T="el'//trim(zoneID)//'", I=',this%solution%interp%M+1, & ', J=',this%solution%interp%M+1 do j = 1,this%solution%interp%M+1 do i = 1,this%solution%interp%M+1 write(fUnit,fmat) x%interior(i,j,iEl,1,1), & x%interior(i,j,iEl,1,2), & solution%interior(i,j,iEl,1:this%solution%nvar), & solutionGradient%interior(i,j,iEl,1:this%solution%nvar,1), & solutionGradient%interior(i,j,iEl,1:this%solution%nvar,2), & dsdt%interior(i,j,iEl,1:this%solution%nvar) enddo enddo enddo close(UNIT=fUnit) call x%Free() call solution%Free() call dsdt%Free() call interp%Free() endsubroutine WriteTecplot_DGModel2D_t endmodule SELF_DGModel2D_t