Merge pull request #85 from FluidNumerics/patches

Patches

CMakeLists.txt

@@ -115,7 +115,7 @@ endif()
 # Check for openmp support if offload target is provided
 if( SELF_ENABLE_MULTITHREADING )
     if( "${CMAKE_Fortran_COMPILER_ID}" STREQUAL "GNU" )
-        set( OFFLOAD_FLAGS "-ftree-parallelize-loops=${SELF_MULTITHREADING_NTHREADS} -fopt-info-loop" )
+        set( OFFLOAD_FLAGS "-ftree-parallelize-loops=${SELF_MULTITHREADING_NTHREADS}" ) #-fopt-info-loop"
     elseif( "${CMAKE_Fortran_COMPILER_ID}" STREQUAL "Intel" )
         set( OFFLOAD_FLAGS "-parallel -qopt-report -qopt-report-phase=par" )
     elseif( "${CMAKE_Fortran_COMPILER_ID}" STREQUAL "IntelLLVM" )

examples/linear_euler2d_spherical_soundwave_closeddomain.f90

@@ -30,6 +30,13 @@ program LinearEuler_Example
   use self_LinearEuler2D
 
   implicit none
+  real(prec),parameter :: rho0 = 1.225_prec
+  real(prec),parameter :: rhoprime = 0.01_prec
+  real(prec),parameter :: c = 1.0_prec ! Speed of sound
+  real(prec),parameter :: Lr = 0.06_prec
+  real(prec),parameter :: x0 = 0.5_prec
+  real(prec),parameter :: y0 = 0.5_prec
+
   character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3'
   integer,parameter :: controlDegree = 7
   integer,parameter :: targetDegree = 15
@@ -64,16 +71,12 @@ program LinearEuler_Example
   ! Initialize the model
   call modelobj%Init(mesh,geometry)
   modelobj%prescribed_bcs_enabled = .false. ! Disables prescribed boundary condition block for gpu accelerated implementations
-  ! modelobj%rho0 = ! optional, set the reference density
-  ! modelobj%c = ! optional set the reference sound wave speed
-  ! modelobj%g = ! optional set the gravitational acceleration (y-direction)
+  modelobj%tecplot_enabled = .false. ! Disables tecplot output
+  modelobj%rho0 = rho0 ! optional, set the reference density
+  modelobj%c = c ! optional set the reference sound wave speed
 
   ! Set the initial condition
-  call modelobj%solution%SetEquation(1,'d = 0.0001*exp( -ln(2)*( (x-0.5)^2 + (y-0.5)^2 )/0.0036 )') ! density
-  call modelobj%solution%SetEquation(2,'u = 0') ! u
-  call modelobj%solution%SetEquation(3,'v = 0') ! v
-  call modelobj%solution%SetEquation(4,'p = 0.0001*exp( -ln(2)*( (x-0.5)^2 + (y-0.5)^2 )/0.0036 )') ! pressure
-  call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
+  call modelobj%SphericalSoundWave(rhoprime,Lr,x0,y0)
 
   call modelobj%WriteModel()
   call modelobj%IncrementIOCounter()
@@ -90,8 +93,8 @@ program LinearEuler_Example
 
   ef = modelobj%entropy
 
-  if(ef > e0) then
-    print*,"Error: Final absmax greater than initial absmax! ",e0,ef
+  if(ef /= ef) then
+    print*,"Error: Final entropy is inf or nan",ef
     stop 1
   endif
   ! Clean up

pyself/model2d.py

@@ -3,15 +3,14 @@
 
 
 # Other SELF modules
-import self.geometry as geometry
+import pyself.geometry as geometry
 
 class model:
     def __init__(self):
         self.solution = None
         self.pvdata = None # Pyvista data
         self.varnames = None
         self.varunits = None
-        self.nvar = 0
         self.geom = geometry.semquad()
 
     def load(self, hdf5File):
@@ -26,26 +25,56 @@ def load(self, hdf5File):
 
         if 'controlgrid' in list(f.keys()):
 
-            s = f['controlgrid/solution']
-            self.solution = []
-            i = 0
-            for k in s.keys():
-                if k == 'metadata':
-                    for v in s[f"{k}/name"].keys():
-                        self.solution.append( {"name":s[f"{k}/name/{v}"][()][0].decode('utf-8'), "units":s[f"{k}/units/{v}"][()][0].decode('utf-8'), 'data': None} )
-                    self.nvar = len(self.solution)
+            controlgrid = f['controlgrid']
+            for group_name in controlgrid.keys():
+
+                if( group_name == 'geometry' ):
+                    continue
+
+                group = controlgrid[group_name]
+                # Create a list to hold data for this group
+                setattr(self, group_name, [])
+                group_data = getattr(self, group_name)
+                print(f"Loading {group_name} group")
+
+                # Load metadata information
+                if( 'metadata' in list(group.keys()) ):
+                    for v in group[f"metadata/name"].keys():
+
+                        name = group[f"metadata/name/{v}"].asstr()[()][0]
+                        try:
+                            units = group[f"metadata/units/{v}"].asstr()[()][0]
+                        except:
+                            units = "error"
+
+                        group_data.append({
+                            "name": name,
+                            "units": units,
+                            'data': None
+                        })
                 else:
-                    d = s[k]
-                    N = d.shape[2]
-                    # Find index for this field
-                    i=0
-                    for sol in self.solution:
-                        if sol['name'] == k:
-                            break
-                        else:
-                            i+=1
-
-                    self.solution[i]['data'] = da.from_array(d, chunks=(self.geom.daskChunkSize,N,N))
+                    print(f"Error: /controlgrid/{group_name}/metadata group not found in {hdf5File}.")
+                    return 1
+
+                for k in group.keys():
+                    k_decoded = k.encode('utf-8').decode('utf-8')
+                    if k == 'metadata':
+                        continue
+                    else:
+                        print(f"Loading {k_decoded} field")
+                        # Load the actual data
+                        d = group[k]
+                        N = d.shape[2]
+
+                        # Find index for this field
+                        i = 0
+                        for sol in group_data:
+                            if sol['name'] == k_decoded:
+                                break
+                            else:
+                                i += 1
+
+                        group_data[i]['data'] = da.from_array(d, chunks=(self.geom.daskChunkSize, N, N))
 
             self.generate_pyvista()
 
@@ -97,8 +126,13 @@ def generate_pyvista(self):
 
         # Load fields into pvdata
         k = 0
-        for var in self.solution:
-            self.pvdata.point_data.set_array(var['data'].flatten(),var['name'])
-            k+=1
+        for attr in self.__dict__:
+            if not attr in ['pvdata','varnames','varunits','geom']:
+                controlgroup = getattr(self, attr)
+                #print(f"Loading {attr} into pvdata")
+                for var in controlgroup:
+                 #   print(f"Loading {var['name']} into pvdata")
+                    self.pvdata.point_data.set_array(var['data'].flatten(),var['name'])
+                    k+=1
 
         print(self.pvdata)

setup.py

@@ -11,7 +11,8 @@
     packages=['pyself'],
     install_requires=['h5py>=3.7.0',
                       'dask',
-                      'pyvista'],
+                      'pyvista',
+                      'imageio[ffmpeg]'],
     classifiers=[
         'Development Status :: 1 - Planning',
         'Intended Audience :: Science/Research',

src/SELF_DGModel2D_t.f90

@@ -64,6 +64,7 @@ module SELF_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
 
@@ -149,6 +150,43 @@ subroutine Free_DGModel2D_t(this)
 
   endsubroutine Free_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
@@ -306,7 +344,7 @@ subroutine CalculateEntropy_DGModel2D_t(this)
     do iel = 1,this%geometry%nelem
       do j = 1,this%solution%interp%N+1
         do i = 1,this%solution%interp%N+1
-          jac = this%geometry%J%interior(i,j,iel,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
         enddo

src/SELF_DGModel3D_t.f90

@@ -62,6 +62,7 @@ module SELF_DGModel3D_t
     procedure :: SourceMethod => sourcemethod_DGModel3D_t
     procedure :: SetBoundaryCondition => setboundarycondition_DGModel3D_t
     procedure :: SetGradientBoundaryCondition => setgradientboundarycondition_DGModel3D_t
+    procedure :: ReportMetrics => ReportMetrics_DGModel3D_t
 
     procedure :: UpdateSolution => UpdateSolution_DGModel3D_t
 
@@ -147,6 +148,43 @@ subroutine Free_DGModel3D_t(this)
 
   endsubroutine Free_DGModel3D_t
 
+  subroutine ReportMetrics_DGModel3D_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(DGModel3D_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_DGModel3D_t
+
   subroutine SetSolutionFromEqn_DGModel3D_t(this,eqn)
     implicit none
     class(DGModel3D_t),intent(inout) :: this
@@ -305,7 +343,7 @@ subroutine CalculateEntropy_DGModel3D_t(this)
       do k = 1,this%solution%interp%N+1
         do j = 1,this%solution%interp%N+1
           do i = 1,this%solution%interp%N+1
-            jac = this%geometry%J%interior(i,j,k,iel,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
           enddo

src/SELF_LinearEuler2D_t.f90

@@ -74,6 +74,7 @@ module self_LinearEuler2D_t
     procedure :: flux2d => flux2d_LinearEuler2D_t
     procedure :: riemannflux2d => riemannflux2d_LinearEuler2D_t
     !procedure :: source2d => source2d_LinearEuler2D_t
+    procedure :: SphericalSoundWave => SphericalSoundWave_LinearEuler2D_t
 
   endtype LinearEuler2D_t
 
@@ -100,7 +101,7 @@ subroutine SetMetadata_LinearEuler2D_t(this)
     call this%solution%SetName(3,"v") ! y-velocity component
     call this%solution%SetUnits(3,"m⋅s⁻¹")
 
-    call this%solution%SetName(4,"pressure") ! Pressure
+    call this%solution%SetName(4,"P") ! Pressure
     call this%solution%SetUnits(4,"kg⋅m⁻¹⋅s⁻²")
 
   endsubroutine SetMetadata_LinearEuler2D_t
@@ -188,4 +189,49 @@ pure function riemannflux2d_LinearEuler2D_t(this,sL,sR,dsdx,nhat) result(flux)
 
   endfunction riemannflux2d_LinearEuler2D_t
 
+  subroutine SphericalSoundWave_LinearEuler2D_t(this,rhoprime,Lr,x0,y0)
+    !! This subroutine sets the initial condition for a weak blast wave
+    !! problem. The initial condition is given by
+    !!
+    !! \begin{equation}
+    !! \begin{aligned}
+    !! \rho &= \rho_0 + \rho' \exp\left( -\ln(2) \frac{(x-x_0)^2 + (y-y_0)^2}{L_r^2} \right)
+    !! u &= 0 \\
+    !! v &= 0 \\
+    !! E &= \frac{P_0}{\gamma - 1} + E \exp\left( -\ln(2) \frac{(x-x_0)^2 + (y-y_0)^2}{L_e^2} \right)
+    !! \end{aligned}
+    !! \end{equation}
+    !!
+    implicit none
+    class(LinearEuler2D_t),intent(inout) :: this
+    real(prec),intent(in) ::  rhoprime,Lr,x0,y0
+    ! Local
+    integer :: i,j,iEl
+    real(prec) :: x,y,rho,r,E
+
+    print*,__FILE__," : Configuring weak blast wave initial condition. "
+    print*,__FILE__," : rhoprime = ",rhoprime
+    print*,__FILE__," : Lr = ",Lr
+    print*,__FILE__," : x0 = ",x0
+    print*,__FILE__," : y0 = ",y0
+
+    do concurrent(i=1:this%solution%N+1,j=1:this%solution%N+1, &
+                  iel=1:this%mesh%nElem)
+      x = this%geometry%x%interior(i,j,iEl,1,1)-x0
+      y = this%geometry%x%interior(i,j,iEl,1,2)-y0
+      r = sqrt(x**2+y**2)
+
+      rho = (rhoprime)*exp(-log(2.0_prec)*r**2/Lr**2)
+
+      this%solution%interior(i,j,iEl,1) = rho
+      this%solution%interior(i,j,iEl,2) = 0.0_prec
+      this%solution%interior(i,j,iEl,3) = 0.0_prec
+      this%solution%interior(i,j,iEl,4) = rho*this%c*this%c
+
+    enddo
+
+    call this%ReportMetrics()
+
+  endsubroutine SphericalSoundWave_LinearEuler2D_t
+
 endmodule self_LinearEuler2D_t

src/SELF_Mesh_2D_t.f90

@@ -225,6 +225,8 @@ subroutine ResetBoundaryConditionType_Mesh2D_t(this,bcid)
       enddo
     enddo
 
+    call this%UpdateDevice()
+
   endsubroutine ResetBoundaryConditionType_Mesh2D_t
 
   subroutine UniformStructuredMesh_Mesh2D_t(this,nxPerTile,nyPerTile,nTileX,nTileY,dx,dy,bcids)

src/SELF_Mesh_3D_t.f90

@@ -301,6 +301,8 @@ subroutine ResetBoundaryConditionType_Mesh3D_t(this,bcid)
       enddo
     enddo
 
+    call this%UpdateDevice()
+
   endsubroutine ResetBoundaryConditionType_Mesh3D_t
 
   subroutine RecalculateFlip_Mesh3D_t(this)