! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! ! ! 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_MultiSetModel3D_t !! A parent model that runs one DG model per element set and couples the sets at their !! shared faces. !! !! The global mesh is split into nSets element sets (SELF_ElementSetMesh_3D). Set k gets !! its own sub-mesh, geometry and DG model, which may be any extension of DGModel3D: for !! example a fluid model on one set and a solid model on another. The parent extends !! Model, so the stock time integrators drive it: every stage calls the parent's !! CalculateTendency and UpdateGRKn, and the parent forwards them to each set. !! !! Setup takes three steps, because user model types cannot be constructed generically: !! !! 1. call ms%Init(globalMesh,interp,setOfElemLocal,nSets) !! builds the sub-meshes and geometries. !! 2. On every rank and for every set k, including sets with no elements on that rank: !! allocate(<ModelType> :: ms%sets(k)%model) !! call ms%sets(k)%model%Init(ms%sets(k)%mesh,ms%sets(k)%geometry) !! then set the initial condition. Model Init is collective over the set's !! communicator (MapBoundaryConditions reduces over it), so no rank may skip a set. !! 3. call ms%SetupInterfaces() !! checks the set models, builds the interface map and registers the interface !! boundary condition on each set model. !! !! Interface faces. In each sub-mesh, a face shared with another set carries !! sideInfo(5) = SELF_BC_SET_INTERFACE. SetupInterfaces registers a no-op hyperbolic !! boundary condition for that id on every set model and maps boundary conditions again, !! so the faces count as mapped and nothing but the parent writes their exterior state. !! !! Each stage, CalculateTendency !! 0. calls the parent's PreTendencyHook, !! 1. interpolates the solution of every set with interface faces on this rank to its !! element faces (BoundaryInterp), !! 2. exchanges the interface traces (SetInterface3D%Exchange), so iface%ghost holds !! the partner set's trace on every local interface face, !! 3. calls FillInterfaceBoundary(k) for each set k, which writes solution%extBoundary !! on set k's interface faces, !! 4. calls each set model's CalculateTendency. !! The default FillInterfaceBoundary copies the ghost trace unchanged and requires both !! sets to have the same nvar. Override it to couple sets whose variables differ. !! !! Invariants: !! - Sets are processed in the same fixed order (k = 1..nSets) on every rank. !! - A set with no elements on a rank skips CalculateTendency and the solution updates !! on that rank. Set models must therefore not call collectives, or post messages on !! their mesh communicator, inside CalculateTendency or anything it calls !! (PreTendencyHook, boundary condition handlers). Every other forwarded procedure !! runs on every set, empty or not. !! - sets(:) must never be reassigned, reallocated or moved: each model stores pointers !! to its own mesh and geometry. !! - Solution gradients are not supported (gradient_enabled must be false on every set): !! the gradient exchange inside each set's CalculateTendency has no hook for interface !! faces. !! - Each sub-mesh owns a duplicate of the global mesh communicator, which Free releases. !! !! Output. WriteModel writes one HDF5 file: !! /interp/... the shared interpolant !! /elementsets/setOfElem set id of every global element !! /set_<k>/controlgrid/solution/<var> (Np,Np,Np,nGlobalElem_k) per variable !! /set_<k>/controlgrid/geometry/<x>_dim<d> node coordinates, d = 1..3 !! with <k> written as I4.4 and elements in set-global order. ReadModel reads the !! solutions back into a model built on the same partition. use SELF_Constants use SELF_Lagrange use SELF_Mesh_3D use SELF_Geometry_3D use SELF_Model use SELF_DGModel3D use SELF_BoundaryConditions use SELF_ElementSetMesh_3D use SELF_SetInterface_3D use SELF_HDF5 use SELF_MultiSetIO_3D use HDF5 use mpi implicit none type ElementSet3D !! One element set: its sub-mesh, geometry and model. The mesh and geometry live on !! the heap so the pointers the model stores to them stay valid. class(DGModel3D),allocatable :: model type(Mesh3D),pointer :: mesh => null() type(SEMHex),pointer :: geometry => null() !! Duplicate of the global mesh communicator, used by this set's decomposition. !! Owned by the parent and freed after the sub-mesh. integer :: comm = MPI_COMM_NULL endtype ElementSet3D type,extends(Model) :: MultiSetModel3D_t integer :: nSets = 0 type(ElementSet3D),allocatable :: sets(:) !! The global mesh the sets were extracted from. Not owned. type(Mesh3D),pointer :: globalMesh => null() !! Set id of every global element, s(e), e = 1..nGlobalElem (replicated). integer,allocatable :: setOfElem(:) !! Set-global id g_k(e) of every global element within its own set (replicated). integer,allocatable :: setGlobalId(:) !! Number of set-k elements on lower ranks, k = 1..nSets. integer,allocatable :: setOffsetThisRank(:) !! Number of interface faces summed over all sets and all ranks. integer :: nInterfaceFaces = 0 !! Interface map and ghost traces; built by SetupInterfaces. type(SetInterface3D) :: iface !! Each set's solution%boundary, the trace the interface exchange sends. type(SetTraceRef),allocatable :: traces(:) logical :: interfacesReady = .false. contains procedure :: Init => Init_MultiSetModel3D_t procedure :: SetupInterfaces => SetupInterfaces_MultiSetModel3D_t procedure :: FillInterfaceBoundary => FillInterfaceBoundary_MultiSetModel3D_t procedure :: Free => Free_MultiSetModel3D_t procedure :: CalculateTendency => CalculateTendency_MultiSetModel3D_t procedure :: UpdateSolution => UpdateSolution_MultiSetModel3D_t procedure :: UpdateGRK2 => UpdateGRK2_MultiSetModel3D_t procedure :: UpdateGRK3 => UpdateGRK3_MultiSetModel3D_t procedure :: UpdateGRK4 => UpdateGRK4_MultiSetModel3D_t procedure :: PreStepHook => PreStepHook_MultiSetModel3D_t procedure :: PostStepHook => PostStepHook_MultiSetModel3D_t procedure :: IncrementIOCounter => IncrementIOCounter_MultiSetModel3D_t procedure :: CalculateEntropy => CalculateEntropy_MultiSetModel3D_t procedure :: ReportMetrics => ReportMetrics_MultiSetModel3D_t procedure :: ReportUnmappedBoundaries => ReportUnmappedBoundaries_MultiSetModel3D_t procedure :: ReadModel => Read_MultiSetModel3D_t procedure :: WriteModel => Write_MultiSetModel3D_t procedure :: WriteTecplot => WriteTecplot_MultiSetModel3D_t procedure,private :: SyncSetClocks endtype MultiSetModel3D_t contains subroutine Init_MultiSetModel3D_t(this,globalMesh,interp,setOfElemLocal,nSets) !! Split globalMesh into nSets element sets and build each set's sub-mesh and !! geometry. The set models are not created here; see the module documentation. !! !! Collective over globalMesh%decomp%mpiComm: every rank must call it. !! !! Input !! - globalMesh : the decomposed global mesh; must outlive this model !! - interp : the interpolant shared by every set's geometry !! - setOfElemLocal(1:globalMesh%nElem) : set id, in 1..nSets, of each rank-local !! element !! - nSets : number of element sets, >= 1 implicit none class(MultiSetModel3D_t),intent(out) :: this type(Mesh3D),intent(in),target :: globalMesh type(Lagrange),pointer,intent(in) :: interp integer,intent(in) :: setOfElemLocal(1:globalMesh%nElem) integer,intent(in) :: nSets ! Local integer :: k,nGlobalElem,ierror integer,allocatable :: setOffsets(:,:) this%globalMesh => globalMesh this%nSets = nSets this%nvar = 0 this%t = 0.0_prec this%dt = 0.0_prec nGlobalElem = globalMesh%decomp%offsetElem(globalMesh%decomp%nRanks+1) allocate(this%setOfElem(1:nGlobalElem),this%setGlobalId(1:nGlobalElem)) call GatherSetOfElem(globalMesh,setOfElemLocal,this%setOfElem) if(.not. ValidateSetPartition(globalMesh,this%setOfElem,nSets)) then error stop "MultiSetModel3D : invalid element-set partition" endif allocate(setOffsets(1:globalMesh%decomp%nRanks+1,1:nSets)) call BuildSetNumbering(globalMesh,this%setOfElem,nSets,this%setGlobalId,setOffsets) this%setOffsetThisRank = setOffsets(globalMesh%decomp%rankId+1,1:nSets) allocate(this%sets(1:nSets)) do k = 1,nSets call mpi_comm_dup(globalMesh%decomp%mpiComm,this%sets(k)%comm,ierror) allocate(this%sets(k)%mesh) call ExtractSubMesh(globalMesh,k,this%setOfElem,this%setGlobalId,setOffsets, & this%sets(k)%comm,this%sets(k)%mesh) allocate(this%sets(k)%geometry) call this%sets(k)%geometry%Init(interp,this%sets(k)%mesh%nElem) ! Element-local, and a zero-element launch is invalid on the GPU, so a set with ! no elements on this rank skips it. if(this%sets(k)%mesh%nElem > 0) then call this%sets(k)%geometry%GenerateFromMesh(this%sets(k)%mesh) endif enddo deallocate(setOffsets) endsubroutine Init_MultiSetModel3D_t subroutine SetupInterfaces_MultiSetModel3D_t(this) !! Finish setup once every set model has been allocated and initialized. !! !! Checks that every set model exists, was initialized on its own set's mesh and !! geometry, and has gradients disabled. Builds the interface map, then registers a no-op hyperbolic !! boundary condition for SELF_BC_SET_INTERFACE on every set model and maps its !! boundary conditions again, so interface faces are not counted as unmapped !! boundaries. !! !! Collective: MapBoundaryConditions reduces over each set's communicator, and the !! interface map and face count use the global one. Every rank must call it. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local procedure(SELF_bcMethod),pointer :: bcfunc integer :: k,iEl,s,nLocal,ierror integer,allocatable :: nvarOfSet(:) do k = 1,this%nSets if(.not. allocated(this%sets(k)%model)) then print*,__FILE__," : set ",k," has no model allocated" error stop "MultiSetModel3D : every set needs a model before SetupInterfaces" endif if(.not. associated(this%sets(k)%model%mesh,this%sets(k)%mesh)) then print*,__FILE__," : the model of set ",k," was not initialized on that set's mesh" error stop "MultiSetModel3D : initialize sets(k)%model on sets(k)%mesh" endif if(.not. associated(this%sets(k)%model%geometry,this%sets(k)%geometry)) then print*,__FILE__," : the model of set ",k," was not initialized on that set's geometry" error stop "MultiSetModel3D : initialize sets(k)%model on sets(k)%geometry" endif if(this%sets(k)%model%gradient_enabled) then print*,__FILE__," : the model of set ",k," has gradient_enabled" error stop "MultiSetModel3D : solution gradients are not supported" endif enddo nLocal = 0 do k = 1,this%nSets do iEl = 1,this%sets(k)%mesh%nElem do s = 1,6 if(this%sets(k)%mesh%sideInfo(3,s,iEl) == 0 .and. & this%sets(k)%mesh%sideInfo(5,s,iEl) == SELF_BC_SET_INTERFACE) then nLocal = nLocal+1 endif enddo enddo enddo call mpi_allreduce(nLocal,this%nInterfaceFaces,1,MPI_INTEGER,MPI_SUM, & this%globalMesh%decomp%mpiComm,ierror) allocate(nvarOfSet(1:this%nSets)) do k = 1,this%nSets nvarOfSet(k) = this%sets(k)%model%nvar enddo call this%iface%BuildMap(this%globalMesh,this%setOfElem,this%setGlobalId, & this%setOffsetThisRank,nvarOfSet, & this%sets(1)%model%solution%interp%N) deallocate(nvarOfSet) ! The sub-meshes tag a face as an interface exactly when the global mesh puts the ! two elements in different sets, so both counts must agree. if(this%iface%nFaces /= nLocal) then error stop "MultiSetModel3D : interface map disagrees with the sub-mesh tags" endif allocate(this%traces(1:this%nSets)) do k = 1,this%nSets this%traces(k)%trace => this%sets(k)%model%solution%boundary enddo ! Size each sub-mesh's request scratch for its own halo before stepping: a send and a ! receive per rank-remote side (or mortar sub-side) per variable, as the exchanges ! reserve. The sub-meshes are built with a minimal scratch. do k = 1,this%nSets call ReserveSetMessages(this%sets(k)%mesh,this%sets(k)%model%nvar) enddo bcfunc => hbc3d_SetInterface_MultiSetModel3D do k = 1,this%nSets call this%sets(k)%model%hyperbolicBCs%RegisterBoundaryCondition( & SELF_BC_SET_INTERFACE,"set_interface",bcfunc) call this%sets(k)%model%MapBoundaryConditions() enddo this%interfacesReady = .true. endsubroutine SetupInterfaces_MultiSetModel3D_t subroutine ReserveSetMessages(mesh,nvar) !! Reserve mesh%decomp's request scratch for the largest scalar exchange a set with !! nvar variables posts: 2*nvar messages per rank-remote conforming side, mortar !! sub-side or anisotropic mortar half. Rank-local; no communication. implicit none type(Mesh3D),intent(inout) :: mesh integer,intent(in) :: nvar ! Local integer :: nRemote nRemote = mesh%decomp%CountRemoteSides(mesh%sideInfo,mesh%nElem,6) if(mesh%nMortars > 0) then nRemote = max(nRemote,mesh%decomp%CountRemoteMortarSides(mesh%mortarInfo,mesh%nMortars,4)) endif if(mesh%nAnisoMortars > 0) then nRemote = max(nRemote, & mesh%decomp%CountRemoteMortarSides(mesh%anisoMortarInfo,mesh%nAnisoMortars,2)) endif call mesh%decomp%ReserveMessages(2*nvar*nRemote) endsubroutine ReserveSetMessages subroutine hbc3d_SetInterface_MultiSetModel3D(bc,mymodel) !! Boundary condition for SELF_BC_SET_INTERFACE faces. It writes nothing: the parent !! fills solution%extBoundary on these faces before each set's CalculateTendency. !! Registering it keeps the faces out of the unmapped-boundary count. implicit none class(BoundaryCondition),intent(in) :: bc class(Model),intent(inout) :: mymodel if(.false.) mymodel%nvar = bc%nBoundaries ! suppress unused-dummy-argument warning endsubroutine hbc3d_SetInterface_MultiSetModel3D subroutine FillInterfaceBoundary_MultiSetModel3D_t(this,iset) !! Write solution%extBoundary on the interface faces of set iset from the ghost !! traces of the last exchange: !! !! extBoundary(i,j,s,e,v) = ghost(i,j,v,f), v = 1..nvar, !! !! for each interface face f = (e,s) of the set. The partner set must have the same !! nvar; otherwise this stops with an error. Override it to convert between the !! variables of different set models (for example acoustic and elastic states). !! Runs on the host; a GPU override must write the device copy of extBoundary. implicit none class(MultiSetModel3D_t),intent(inout) :: this integer,intent(in) :: iset ! Local integer :: f,nvar nvar = this%sets(iset)%model%nvar do f = this%iface%setFaceStart(iset),this%iface%setFaceStart(iset+1)-1 if(this%iface%nvarOfSet(this%iface%facePartnerSet(f)) /= nvar) then print*,__FILE__," : sets ",iset," and ",this%iface%facePartnerSet(f), & " have different nvar" error stop "MultiSetModel3D : override FillInterfaceBoundary to couple these sets" endif this%sets(iset)%model%solution%extBoundary(:,:,this%iface%faceSide(f), & this%iface%faceElem(f),1:nvar) = & this%iface%ghost(:,:,1:nvar,f) enddo endsubroutine FillInterfaceBoundary_MultiSetModel3D_t subroutine Free_MultiSetModel3D_t(this) !! Free the interface map, then every set's model, geometry and sub-mesh, then the !! sub-mesh communicators. The global mesh is not owned and is left alone. !! Collective. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k,ierror call this%iface%Free() if(allocated(this%traces)) deallocate(this%traces) if(allocated(this%sets)) then do k = 1,this%nSets if(allocated(this%sets(k)%model)) then call this%sets(k)%model%Free() deallocate(this%sets(k)%model) endif if(associated(this%sets(k)%geometry)) then call this%sets(k)%geometry%Free() deallocate(this%sets(k)%geometry) endif if(associated(this%sets(k)%mesh)) then call this%sets(k)%mesh%Free() deallocate(this%sets(k)%mesh) endif if(this%sets(k)%comm /= MPI_COMM_NULL) then call mpi_comm_free(this%sets(k)%comm,ierror) endif enddo deallocate(this%sets) endif if(allocated(this%setOfElem)) deallocate(this%setOfElem) if(allocated(this%setGlobalId)) deallocate(this%setGlobalId) if(allocated(this%setOffsetThisRank)) deallocate(this%setOffsetThisRank) nullify(this%globalMesh) this%nSets = 0 this%nInterfaceFaces = 0 this%interfacesReady = .false. endsubroutine Free_MultiSetModel3D_t subroutine SyncSetClocks(this) !! Copy the parent's simulation time and time step to every set model. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k do k = 1,this%nSets this%sets(k)%model%t = this%t this%sets(k)%model%dt = this%dt enddo endsubroutine SyncSetClocks subroutine CalculateTendency_MultiSetModel3D_t(this) !! Exchange interface traces, fill the interface exterior states, then compute the !! tendency of every set with elements on this rank, in set order. See the module !! documentation for the stage order. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k if(.not. this%interfacesReady) then error stop "MultiSetModel3D : call SetupInterfaces before stepping" endif call this%SyncSetClocks() call this%PreTendencyHook() ! User-supplied; runs before the interface traces are gathered if(this%nInterfaceFaces > 0) then do k = 1,this%nSets if(this%iface%setFaceStart(k+1) > this%iface%setFaceStart(k)) then call this%sets(k)%model%solution%BoundaryInterp() endif enddo call this%iface%Exchange(this%traces) do k = 1,this%nSets if(this%iface%setFaceStart(k+1) > this%iface%setFaceStart(k)) then call this%FillInterfaceBoundary(k) endif enddo endif do k = 1,this%nSets if(this%sets(k)%mesh%nElem > 0) call this%sets(k)%model%CalculateTendency() enddo endsubroutine CalculateTendency_MultiSetModel3D_t subroutine UpdateSolution_MultiSetModel3D_t(this,dt) !! Forward-Euler update of every set with elements on this rank. dt defaults to the !! parent's time step. implicit none class(MultiSetModel3D_t),intent(inout) :: this real(prec),optional,intent(in) :: dt ! Local integer :: k call this%SyncSetClocks() do k = 1,this%nSets if(this%sets(k)%mesh%nElem == 0) cycle if(present(dt)) then call this%sets(k)%model%UpdateSolution(dt) else call this%sets(k)%model%UpdateSolution() endif enddo endsubroutine UpdateSolution_MultiSetModel3D_t subroutine UpdateGRK2_MultiSetModel3D_t(this,m) !! Stage m of the low-storage RK2 update, for every set with elements on this rank. implicit none class(MultiSetModel3D_t),intent(inout) :: this integer,intent(in) :: m ! Local integer :: k call this%SyncSetClocks() do k = 1,this%nSets if(this%sets(k)%mesh%nElem > 0) call this%sets(k)%model%UpdateGRK2(m) enddo endsubroutine UpdateGRK2_MultiSetModel3D_t subroutine UpdateGRK3_MultiSetModel3D_t(this,m) !! Stage m of the low-storage RK3 update, for every set with elements on this rank. implicit none class(MultiSetModel3D_t),intent(inout) :: this integer,intent(in) :: m ! Local integer :: k call this%SyncSetClocks() do k = 1,this%nSets if(this%sets(k)%mesh%nElem > 0) call this%sets(k)%model%UpdateGRK3(m) enddo endsubroutine UpdateGRK3_MultiSetModel3D_t subroutine UpdateGRK4_MultiSetModel3D_t(this,m) !! Stage m of the low-storage RK4 update, for every set with elements on this rank. implicit none class(MultiSetModel3D_t),intent(inout) :: this integer,intent(in) :: m ! Local integer :: k call this%SyncSetClocks() do k = 1,this%nSets if(this%sets(k)%mesh%nElem > 0) call this%sets(k)%model%UpdateGRK4(m) enddo endsubroutine UpdateGRK4_MultiSetModel3D_t subroutine PreStepHook_MultiSetModel3D_t(this) !! Run every set's PreStepHook, empty sets included, after syncing the set clocks. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k call this%SyncSetClocks() do k = 1,this%nSets call this%sets(k)%model%PreStepHook() enddo endsubroutine PreStepHook_MultiSetModel3D_t subroutine PostStepHook_MultiSetModel3D_t(this) !! Run every set's PostStepHook, empty sets included, after syncing the set clocks. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k call this%SyncSetClocks() do k = 1,this%nSets call this%sets(k)%model%PostStepHook() enddo endsubroutine PostStepHook_MultiSetModel3D_t subroutine IncrementIOCounter_MultiSetModel3D_t(this) !! Advance the output counter of the parent and of every set. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k this%ioIterate = this%ioIterate+1 do k = 1,this%nSets call this%sets(k)%model%IncrementIOCounter() enddo endsubroutine IncrementIOCounter_MultiSetModel3D_t subroutine CalculateEntropy_MultiSetModel3D_t(this) !! entropy = sum over sets k of the entropy of set k. Each set's CalculateEntropy is !! collective over its communicator, so it runs on every set, empty or not. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k this%entropy = 0.0_prec do k = 1,this%nSets call this%sets(k)%model%CalculateEntropy() this%entropy = this%entropy+this%sets(k)%model%entropy enddo endsubroutine CalculateEntropy_MultiSetModel3D_t subroutine ReportMetrics_MultiSetModel3D_t(this) !! Report each set's metrics. Rank-local; sets with no elements on this rank are !! skipped, since their min and max are undefined. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k do k = 1,this%nSets if(this%sets(k)%mesh%nElem > 0) call this%sets(k)%model%ReportMetrics() enddo endsubroutine ReportMetrics_MultiSetModel3D_t subroutine ReportUnmappedBoundaries_MultiSetModel3D_t(this) !! Report unmapped boundary faces for every set. The counts are already global. implicit none class(MultiSetModel3D_t),intent(inout) :: this ! Local integer :: k do k = 1,this%nSets call this%sets(k)%model%ReportUnmappedBoundaries() enddo endsubroutine ReportUnmappedBoundaries_MultiSetModel3D_t subroutine Write_MultiSetModel3D_t(this,filename) !! Write every set's solution and geometry to one HDF5 file; see the module !! documentation for the layout. Without filename the file is !! solution.<ioIterate>.h5, as for a single model. !! !! Collective over the global mesh communicator. Ranks with no elements of a set !! still take part in every dataset of that set. implicit none class(MultiSetModel3D_t),intent(inout) :: this character(*),intent(in),optional :: filename ! Local integer(HID_T) :: fileId character(LEN=self_FileNameLength) :: pickupFile character(13) :: timeStampString character(4) :: setString character(1) :: dimString character(LEN=:),allocatable :: group integer :: k,ivar,idir,elemOffset logical :: parallel if(present(filename)) then pickupFile = filename else write(timeStampString,'(I13.13)') this%ioIterate pickupFile = 'solution.'//timeStampString//'.h5' endif parallel = this%globalMesh%decomp%mpiEnabled if(this%globalMesh%decomp%rankId == 0) then print*,__FILE__//" : Writing pickup file : "//trim(pickupFile) endif if(parallel) then call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId,this%globalMesh%decomp%mpiComm) else call Open_HDF5(pickupFile,H5F_ACC_TRUNC_F,fileId) endif call this%sets(1)%model%solution%interp%WriteHDF5(fileId) call CreateGroup_HDF5(fileId,'/elementsets') call WriteArray_HDF5(fileId,'/elementsets/setOfElem',this%setOfElem) do k = 1,this%nSets write(setString,'(I4.4)') k group = '/set_'//setString call CreateGroup_HDF5(fileId,group) call CreateGroup_HDF5(fileId,group//'/controlgrid') call CreateGroup_HDF5(fileId,group//'/controlgrid/solution') call CreateGroup_HDF5(fileId,group//'/controlgrid/geometry') elemOffset = this%sets(k)%mesh%decomp%offsetElem(this%globalMesh%decomp%rankId+1) call this%sets(k)%model%solution%UpdateHost() do ivar = 1,this%sets(k)%model%solution%nVar call WriteSetArray_HDF5(fileId, & group//'/controlgrid/solution/'// & trim(this%sets(k)%model%solution%meta(ivar)%name), & this%sets(k)%model%solution%interior(:,:,:,:,ivar), & elemOffset,this%sets(k)%mesh%nGlobalElem,parallel) enddo do idir = 1,3 write(dimString,'(I1)') idir call WriteSetArray_HDF5(fileId, & group//'/controlgrid/geometry/'// & trim(this%sets(k)%geometry%x%meta(1)%name)//'_dim'//dimString, & this%sets(k)%geometry%x%interior(:,:,:,:,1,idir), & elemOffset,this%sets(k)%mesh%nGlobalElem,parallel) enddo enddo call Close_HDF5(fileId) endsubroutine Write_MultiSetModel3D_t subroutine Read_MultiSetModel3D_t(this,filename) !! Read every set's solution from a file written by WriteModel, then copy it to the !! device. The file must hold the same element-set partition (setOfElem) and the !! same number of nodes per element; a variable with no dataset keeps its current !! value. The host copy is refreshed from the device first, so on GPU builds that !! value is the device's. !! !! Collective over the global mesh communicator. implicit none class(MultiSetModel3D_t),intent(inout) :: this character(*),intent(in) :: filename ! Local integer(HID_T) :: fileId character(4) :: setString character(LEN=:),allocatable :: dsetName integer,allocatable :: fileSetOfElem(:) integer :: k,ivar,elemOffset logical :: parallel parallel = this%globalMesh%decomp%mpiEnabled if(parallel) then call Open_HDF5(filename,H5F_ACC_RDONLY_F,fileId,this%globalMesh%decomp%mpiComm) else call Open_HDF5(filename,H5F_ACC_RDONLY_F,fileId) endif ! Check the stored extent first: ReadArray_HDF5 reads the whole dataset into the ! buffer, so a longer dataset would write past it. if(DatasetLength1D_HDF5(fileId,'/elementsets/setOfElem') /= size(this%setOfElem)) then print*,__FILE__," : ",trim(filename)," holds no element-set partition for ", & size(this%setOfElem)," elements" error stop "MultiSetModel3D : ReadModel partition mismatch" endif allocate(fileSetOfElem(1:size(this%setOfElem))) call ReadArray_HDF5(fileId,'/elementsets/setOfElem',fileSetOfElem) if(any(fileSetOfElem /= this%setOfElem)) then print*,__FILE__," : ",trim(filename)," holds a different element-set partition" error stop "MultiSetModel3D : ReadModel partition mismatch" endif deallocate(fileSetOfElem) do k = 1,this%nSets call this%sets(k)%model%solution%UpdateHost() write(setString,'(I4.4)') k elemOffset = this%sets(k)%mesh%decomp%offsetElem(this%globalMesh%decomp%rankId+1) do ivar = 1,this%sets(k)%model%solution%nVar dsetName = '/set_'//setString//'/controlgrid/solution/'// & trim(this%sets(k)%model%solution%meta(ivar)%name) if(.not. DatasetExists_HDF5(fileId,dsetName)) then if(this%globalMesh%decomp%rankId == 0) then print*,__FILE__," : Pickup file holds no ",dsetName, & " - keeping the current value for this variable." endif cycle endif call ReadSetArray_HDF5(fileId,dsetName, & this%sets(k)%model%solution%interior(:,:,:,:,ivar), & elemOffset,this%sets(k)%mesh%nGlobalElem,parallel) enddo call this%sets(k)%model%solution%UpdateDevice() enddo call Close_HDF5(fileId) endsubroutine Read_MultiSetModel3D_t subroutine WriteTecplot_MultiSetModel3D_t(this,filename) !! Write one rank-local Tecplot file per set with elements on this rank. !! !! Without filename, set k writes set<k>.solution.<rank>.<ioIterate>.tec. With !! filename, set k writes set<k>.<filename>. Rank-local; no collectives. implicit none class(MultiSetModel3D_t),intent(inout) :: this character(*),intent(in),optional :: filename ! Local integer :: k character(5) :: rankString character(13) :: timeStampString character(4) :: setString character(LEN=self_FileNameLength) :: tecFile write(rankString,'(I5.5)') this%globalMesh%decomp%rankId write(timeStampString,'(I13.13)') this%ioIterate do k = 1,this%nSets if(this%sets(k)%mesh%nElem == 0) cycle write(setString,'(I4.4)') k if(present(filename)) then tecFile = 'set'//setString//'.'//trim(filename) else tecFile = 'set'//setString//'.solution.'//rankString//'.'//timeStampString//'.tec' endif call this%sets(k)%model%WriteTecplot(trim(tecFile)) enddo endsubroutine WriteTecplot_MultiSetModel3D_t endmodule SELF_MultiSetModel3D_t