Add user metrics reporting

Author: fluidnumerics-joe

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,42 @@ 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
+    end do
+
+    call this%ReportUserMetrics()
+
+  end subroutine ReportMetrics_DGModel2D_t
+
   subroutine SetSolutionFromEqn_DGModel2D_t(this,eqn)
     implicit none
     class(DGModel2D_t),intent(inout) :: this

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,42 @@ 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
+    end do
+
+    call this%ReportUserMetrics()
+
+  end subroutine ReportMetrics_DGModel3D_t
+
   subroutine SetSolutionFromEqn_DGModel3D_t(this,eqn)
     implicit none
     class(DGModel3D_t),intent(inout) :: this

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
 
@@ -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
+
+      end do
+
+      call this%ReportMetrics()
+
+   end subroutine SphericalSoundWave_LinearEuler2D_t
+
 endmodule self_LinearEuler2D_t