From bb42d746875a6069a7381b12e9b1846d76402e15 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 01:21:39 +0000 Subject: [PATCH 01/11] Add user metrics reporting --- src/SELF_DGModel2D_t.f90 | 37 + src/SELF_DGModel3D_t.f90 | 37 + src/SELF_LinearEuler2D_t.f90 | 46 ++ src/SELF_Model.f90 | 1516 ++++++++++++++++++---------------- 4 files changed, 902 insertions(+), 734 deletions(-) diff --git a/src/SELF_DGModel2D_t.f90 b/src/SELF_DGModel2D_t.f90 index 11225f0ae..394262234 100644 --- a/src/SELF_DGModel2D_t.f90 +++ b/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 diff --git a/src/SELF_DGModel3D_t.f90 b/src/SELF_DGModel3D_t.f90 index cefd73e3a..7d818c1f4 100644 --- a/src/SELF_DGModel3D_t.f90 +++ b/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 diff --git a/src/SELF_LinearEuler2D_t.f90 b/src/SELF_LinearEuler2D_t.f90 index 75fdd853b..c5cfdcf1f 100644 --- a/src/SELF_LinearEuler2D_t.f90 +++ b/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 diff --git a/src/SELF_Model.f90 b/src/SELF_Model.f90 index 11e6e23d1..8fe276bf1 100644 --- a/src/SELF_Model.f90 +++ b/src/SELF_Model.f90 @@ -26,257 +26,283 @@ module SELF_Model - use SELF_SupportRoutines - use SELF_Metadata - use SELF_HDF5 - use HDF5 - use FEQParse + use SELF_SupportRoutines + use SELF_Metadata + use SELF_HDF5 + use HDF5 + use FEQParse - implicit none + implicit none #include "SELF_Macros.h" ! //////////////////////////////////////////////// ! ! Time integration parameters - ! Runge-Kutta 2nd Order (Low Storage) - real(prec),parameter :: rk2_a(1:2) = (/0.0_prec,-0.5_prec/) - real(prec),parameter :: rk2_b(1:2) = (/0.5_prec,0.5_prec/) - real(prec),parameter :: rk2_g(1:2) = (/0.5_prec,1.0_prec/) - - ! Williamson's Runge-Kutta 3rd Order (Low Storage) - real(prec),parameter :: rk3_a(1:3) = (/0.0_prec,-5.0_prec/9.0_prec,-153.0_prec/128.0_prec/) - real(prec),parameter :: rk3_b(1:3) = (/0.0_prec,1.0_prec/3.0_prec,3.0_prec/4.0_prec/) - real(prec),parameter :: rk3_g(1:3) = (/1.0_prec/3.0_prec,15.0_prec/16.0_prec,8.0_prec/15.0_prec/) - - ! Carpenter-Kennedy Runge-Kuttta 4th Order (Low Storage) - real(prec),parameter :: rk4_a(1:5) = (/0.0_prec, & - -1.0_prec, & - -1.0_prec/3.0_prec+ & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - -2.0_prec**(1.0_prec/3.0_prec)- & - 2.0_prec**(2.0_prec/3.0_prec)-2.0_prec, & - -1.0_prec+2.0_prec**(1.0_prec/3.0_prec)/) - - real(prec),parameter :: rk4_b(1:5) = (/ & - 0.0_prec, & - 2.0_prec/3.0_prec+2.0_prec**(1.0_prec/3.0_prec)/3.0_prec+ & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - 2.0_prec/3.0_prec+2.0_prec**(1.0_prec/3.0_prec)/3.0_prec+ & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - 1.0_prec/3.0_prec-2.0_prec**(1.0_prec/3.0_prec)/3.0_prec- & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - 1.0_prec/) - - real(prec),parameter :: rk4_g(1:5) = (/ & - 2.0_prec/3.0_prec+2.0_prec**(1.0_prec/3.0_prec)/3.0_prec+ & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - -2.0_prec**(2.0_prec/3.0_prec)/6.0_prec+1.0_prec/6.0_prec, & - -1.0_prec/3.0_prec-2.0_prec*2.0_prec**(1.0_prec/3.0_prec)/3.0_prec- & - 2.0_prec**(2.0_prec/3.0_prec)/3.0_prec, & - 1.0_prec/3.0_prec-2.0_prec**(1.0_prec/3.0_prec)/3.0_prec- & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - 1.0_prec/3.0_prec+2.0_prec**(1.0_prec/3.0_prec)/6.0_prec+ & - 2.0_prec**(2.0_prec/3.0_prec)/12.0_prec/) + ! Runge-Kutta 2nd Order (Low Storage) + real(prec), parameter :: rk2_a(1:2) = (/0.0_prec, -0.5_prec/) + real(prec), parameter :: rk2_b(1:2) = (/0.5_prec, 0.5_prec/) + real(prec), parameter :: rk2_g(1:2) = (/0.5_prec, 1.0_prec/) + + ! Williamson's Runge-Kutta 3rd Order (Low Storage) + real(prec), parameter :: rk3_a(1:3) = (/0.0_prec, -5.0_prec/9.0_prec, -153.0_prec/128.0_prec/) + real(prec), parameter :: rk3_b(1:3) = (/0.0_prec, 1.0_prec/3.0_prec, 3.0_prec/4.0_prec/) + real(prec), parameter :: rk3_g(1:3) = (/1.0_prec/3.0_prec, 15.0_prec/16.0_prec, 8.0_prec/15.0_prec/) + + ! Carpenter-Kennedy Runge-Kuttta 4th Order (Low Storage) + real(prec), parameter :: rk4_a(1:5) = (/0.0_prec, & + -1.0_prec, & + -1.0_prec/3.0_prec + & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + -2.0_prec**(1.0_prec/3.0_prec) - & + 2.0_prec**(2.0_prec/3.0_prec) - 2.0_prec, & + -1.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/) + + real(prec), parameter :: rk4_b(1:5) = (/ & + 0.0_prec, & + 2.0_prec/3.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec + & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + 2.0_prec/3.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec + & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + 1.0_prec/3.0_prec - 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec - & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + 1.0_prec/) + + real(prec), parameter :: rk4_g(1:5) = (/ & + 2.0_prec/3.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec + & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + -2.0_prec**(2.0_prec/3.0_prec)/6.0_prec + 1.0_prec/6.0_prec, & + -1.0_prec/3.0_prec - 2.0_prec*2.0_prec**(1.0_prec/3.0_prec)/3.0_prec - & + 2.0_prec**(2.0_prec/3.0_prec)/3.0_prec, & + 1.0_prec/3.0_prec - 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec - & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + 1.0_prec/3.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/6.0_prec + & + 2.0_prec**(2.0_prec/3.0_prec)/12.0_prec/) ! - integer,parameter :: SELF_EULER = 100 - integer,parameter :: SELF_RK2 = 200 - integer,parameter :: SELF_RK3 = 300 - integer,parameter :: SELF_RK4 = 400 - ! integer,parameter :: SELF_AB2 = 201 - ! integer,parameter :: SELF_AB3 = 301 - ! integer,parameter :: SELF_AB4 = 401 + integer, parameter :: SELF_EULER = 100 + integer, parameter :: SELF_RK2 = 200 + integer, parameter :: SELF_RK3 = 300 + integer, parameter :: SELF_RK4 = 400 + ! integer,parameter :: SELF_AB2 = 201 + ! integer,parameter :: SELF_AB3 = 301 + ! integer,parameter :: SELF_AB4 = 401 - integer,parameter :: SELF_INTEGRATOR_LENGTH = 10 ! max length of integrator methods when specified as char - integer,parameter :: SELF_EQUATION_LENGTH = 500 + integer, parameter :: SELF_INTEGRATOR_LENGTH = 10 ! max length of integrator methods when specified as char + integer, parameter :: SELF_EQUATION_LENGTH = 500 ! //////////////////////////////////////////////// ! ! Model Formulations ! - integer,parameter :: SELF_FORMULATION_LENGTH = 30 ! max length of integrator methods when specified as char + integer, parameter :: SELF_FORMULATION_LENGTH = 30 ! max length of integrator methods when specified as char + + type, abstract :: Model + + ! Time integration attributes + procedure(SELF_timeIntegrator), pointer :: timeIntegrator => Euler_timeIntegrator + real(prec) :: dt + real(prec) :: t + integer :: ioIterate = 0 + logical :: gradient_enabled = .false. + logical :: prescribed_bcs_enabled = .true. + logical :: tecplot_enabled = .true. + integer :: nvar + ! Standard Diagnostics + real(prec) :: entropy ! Mathematical entropy function for the model + + contains + + procedure :: IncrementIOCounter + + procedure :: PrintType => PrintType_Model + + procedure :: SetNumberOfVariables => SetNumberOfVariables_Model + + procedure :: AdditionalInit => AdditionalInit_Model + procedure :: AdditionalFree => AdditionalFree_Model + procedure :: AdditionalOutput => AdditionalOutput_Model + + procedure :: ForwardStep => ForwardStep_Model + + procedure :: Euler_timeIntegrator + + ! Runge-Kutta methods + procedure :: LowStorageRK2_timeIntegrator + procedure(UpdateGRK), deferred :: UpdateGRK2 + + procedure :: LowStorageRK3_timeIntegrator + procedure(UpdateGRK), deferred :: UpdateGRK3 + + procedure :: LowStorageRK4_timeIntegrator + procedure(UpdateGRK), deferred :: UpdateGRK4 + + procedure :: PreTendency => PreTendency_Model + procedure :: entropy_func => entropy_func_Model + + procedure :: flux1D => flux1d_Model + procedure :: flux2D => flux2d_Model + procedure :: flux3D => flux3d_Model + + procedure :: riemannflux1d => riemannflux1d_Model + procedure :: riemannflux2d => riemannflux2d_Model + procedure :: riemannflux3d => riemannflux3d_Model + + procedure :: source1d => source1d_Model + procedure :: source2d => source2d_Model + procedure :: source3d => source3d_Model + + ! Boundary condition functions (hyperbolic) + procedure :: hbc1d_Prescribed => hbc1d_Prescribed_Model + procedure :: hbc1d_Radiation => hbc1d_Generic_Model + procedure :: hbc1d_NoNormalFlow => hbc1d_Generic_Model + procedure :: hbc2d_Prescribed => hbc2d_Prescribed_Model + procedure :: hbc2d_Radiation => hbc2d_Generic_Model + procedure :: hbc2d_NoNormalFlow => hbc2d_Generic_Model + procedure :: hbc3d_Prescribed => hbc3d_Prescribed_Model + procedure :: hbc3d_Radiation => hbc3d_Generic_Model + procedure :: hbc3d_NoNormalFlow => hbc3d_Generic_Model + + ! Boundary condition functions (parabolic) + procedure :: pbc1d_Prescribed => pbc1d_Prescribed_Model + procedure :: pbc1d_Radiation => pbc1d_Generic_Model + procedure :: pbc1d_NoNormalFlow => pbc1d_Generic_Model + procedure :: pbc2d_Prescribed => pbc2d_Prescribed_Model + procedure :: pbc2d_Radiation => pbc2d_Generic_Model + procedure :: pbc2d_NoNormalFlow => pbc2d_Generic_Model + procedure :: pbc3d_Prescribed => pbc3d_Prescribed_Model + procedure :: pbc3d_Radiation => pbc3d_Generic_Model + procedure :: pbc3d_NoNormalFlow => pbc3d_Generic_Model + + procedure :: ReportEntropy => ReportEntropy_Model + procedure :: ReportMetrics => ReportMetrics_Model + procedure :: ReportUserMetrics => ReportUserMetrics_Model + procedure :: CalculateEntropy => CalculateEntropy_Model + + procedure(UpdateSolution), deferred :: UpdateSolution + procedure(CalculateTendency), deferred :: CalculateTendency + procedure(ReadModel), deferred :: ReadModel + procedure(WriteModel), deferred :: WriteModel + procedure(WriteTecplot), deferred :: WriteTecplot + + generic :: SetTimeIntegrator => SetTimeIntegrator_withChar + procedure, private :: SetTimeIntegrator_withChar + + procedure :: SetSimulationTime + procedure :: GetSimulationTime + + end type Model + + interface + subroutine SELF_timeIntegrator(this, tn) + use SELF_Constants, only: prec + import Model + implicit none + class(Model), intent(inout) :: this + real(prec), intent(in) :: tn + end subroutine SELF_timeIntegrator + end interface + + interface + subroutine UpdateGRK(this, m) + import Model + implicit none + class(Model), intent(inout) :: this + integer, intent(in) :: m + end subroutine UpdateGRK + end interface + + interface + subroutine UpdateSolution(this, dt) + use SELF_Constants, only: prec + import Model + implicit none + class(Model), intent(inout) :: this + real(prec), optional, intent(in) :: dt + end subroutine UpdateSolution + end interface + + interface + subroutine CalculateTendency(this) + import Model + implicit none + class(Model), intent(inout) :: this + end subroutine CalculateTendency + end interface + + interface + subroutine WriteModel(this, filename) + import Model + implicit none + class(Model), intent(inout) :: this + character(*), intent(in), optional :: filename + end subroutine WriteModel + end interface + + interface + subroutine ReadModel(this, filename) + import Model + implicit none + class(Model), intent(inout) :: this + character(*), intent(in) :: filename + end subroutine ReadModel + end interface + + interface + subroutine WriteTecplot(this, filename) + import Model + implicit none + class(Model), intent(inout) :: this + character(*), intent(in), optional :: filename + end subroutine WriteTecplot + end interface - type,abstract :: Model - - ! Time integration attributes - procedure(SELF_timeIntegrator),pointer :: timeIntegrator => Euler_timeIntegrator - real(prec) :: dt - real(prec) :: t - integer :: ioIterate = 0 - logical :: gradient_enabled = .false. - logical :: prescribed_bcs_enabled = .true. - integer :: nvar - ! Standard Diagnostics - real(prec) :: entropy ! Mathematical entropy function for the model - - contains - - procedure :: IncrementIOCounter - - procedure :: PrintType => PrintType_Model - - procedure :: SetNumberOfVariables => SetNumberOfVariables_Model - - procedure :: ForwardStep => ForwardStep_Model - - procedure :: Euler_timeIntegrator - - ! Runge-Kutta methods - procedure :: LowStorageRK2_timeIntegrator - procedure(UpdateGRK),deferred :: UpdateGRK2 - - procedure :: LowStorageRK3_timeIntegrator - procedure(UpdateGRK),deferred :: UpdateGRK3 - - procedure :: LowStorageRK4_timeIntegrator - procedure(UpdateGRK),deferred :: UpdateGRK4 - - procedure :: PreTendency => PreTendency_Model - procedure :: entropy_func => entropy_func_Model - - procedure :: flux1D => flux1d_Model - procedure :: flux2D => flux2d_Model - procedure :: flux3D => flux3d_Model - - procedure :: riemannflux1d => riemannflux1d_Model - procedure :: riemannflux2d => riemannflux2d_Model - procedure :: riemannflux3d => riemannflux3d_Model - - procedure :: source1d => source1d_Model - procedure :: source2d => source2d_Model - procedure :: source3d => source3d_Model - - ! Boundary condition functions (hyperbolic) - procedure :: hbc1d_Prescribed => hbc1d_Prescribed_Model - procedure :: hbc1d_Radiation => hbc1d_Generic_Model - procedure :: hbc1d_NoNormalFlow => hbc1d_Generic_Model - procedure :: hbc2d_Prescribed => hbc2d_Prescribed_Model - procedure :: hbc2d_Radiation => hbc2d_Generic_Model - procedure :: hbc2d_NoNormalFlow => hbc2d_Generic_Model - procedure :: hbc3d_Prescribed => hbc3d_Prescribed_Model - procedure :: hbc3d_Radiation => hbc3d_Generic_Model - procedure :: hbc3d_NoNormalFlow => hbc3d_Generic_Model +contains - ! Boundary condition functions (parabolic) - procedure :: pbc1d_Prescribed => pbc1d_Prescribed_Model - procedure :: pbc1d_Radiation => pbc1d_Generic_Model - procedure :: pbc1d_NoNormalFlow => pbc1d_Generic_Model - procedure :: pbc2d_Prescribed => pbc2d_Prescribed_Model - procedure :: pbc2d_Radiation => pbc2d_Generic_Model - procedure :: pbc2d_NoNormalFlow => pbc2d_Generic_Model - procedure :: pbc3d_Prescribed => pbc3d_Prescribed_Model - procedure :: pbc3d_Radiation => pbc3d_Generic_Model - procedure :: pbc3d_NoNormalFlow => pbc3d_Generic_Model + subroutine IncrementIOCounter(this) + implicit none + class(Model), intent(inout) :: this - procedure :: ReportEntropy => ReportEntropy_Model - procedure :: CalculateEntropy => CalculateEntropy_Model + ! Increment the ioIterate + this%ioIterate = this%ioIterate + 1 - procedure(UpdateSolution),deferred :: UpdateSolution - procedure(CalculateTendency),deferred :: CalculateTendency - procedure(ReadModel),deferred :: ReadModel - procedure(WriteModel),deferred :: WriteModel - procedure(WriteTecplot),deferred :: WriteTecplot + end subroutine IncrementIOCounter - generic :: SetTimeIntegrator => SetTimeIntegrator_withChar - procedure,private :: SetTimeIntegrator_withChar + subroutine SetNumberOfVariables_Model(this) + implicit none + class(Model), intent(inout) :: this - procedure :: SetSimulationTime - procedure :: GetSimulationTime + this%nvar = 1 - endtype Model + end subroutine SetNumberOfVariables_Model - interface - subroutine SELF_timeIntegrator(this,tn) - use SELF_Constants,only:prec - import Model - implicit none - class(Model),intent(inout) :: this - real(prec),intent(in) :: tn - endsubroutine SELF_timeIntegrator - endinterface - - interface - subroutine UpdateGRK(this,m) - import Model - implicit none - class(Model),intent(inout) :: this - integer,intent(in) :: m - endsubroutine UpdateGRK - endinterface - - interface - subroutine UpdateSolution(this,dt) - use SELF_Constants,only:prec - import Model + subroutine AdditionalInit_Model(this) implicit none - class(Model),intent(inout) :: this - real(prec),optional,intent(in) :: dt - endsubroutine UpdateSolution - endinterface - - interface - subroutine CalculateTendency(this) - import Model - implicit none - class(Model),intent(inout) :: this - endsubroutine CalculateTendency - endinterface + class(Model), intent(inout) :: this + return + end subroutine AdditionalInit_Model - interface - subroutine WriteModel(this,filename) - import Model - implicit none - class(Model),intent(inout) :: this - character(*),intent(in),optional :: filename - endsubroutine WriteModel - endinterface - - interface - subroutine ReadModel(this,filename) - import Model + subroutine AdditionalFree_Model(this) implicit none - class(Model),intent(inout) :: this - character(*),intent(in) :: filename - endsubroutine ReadModel - endinterface - - interface - subroutine WriteTecplot(this,filename) - import Model - implicit none - class(Model),intent(inout) :: this - character(*),intent(in),optional :: filename - endsubroutine WriteTecplot - endinterface - -contains - - subroutine IncrementIOCounter(this) - implicit none - class(Model),intent(inout) :: this + class(Model), intent(inout) :: this + return + end subroutine AdditionalFree_Model - ! Increment the ioIterate - this%ioIterate = this%ioIterate+1 - - endsubroutine IncrementIOCounter - - subroutine SetNumberOfVariables_Model(this) - implicit none - class(Model),intent(inout) :: this - - this%nvar = 1 - - endsubroutine SetNumberOfVariables_Model + subroutine AdditionalOutput_Model(this, fileid) + implicit none + class(Model), intent(inout) :: this + integer(HID_T), intent(in) :: fileid + return + end subroutine AdditionalOutput_Model - subroutine PrintType_Model(this) - implicit none - class(Model),intent(in) :: this + subroutine PrintType_Model(this) + implicit none + class(Model), intent(in) :: this - print*,__FILE__//" : Model : No model type" + print *, __FILE__//" : Model : No model type" - endsubroutine PrintType_Model + end subroutine PrintType_Model - subroutine PreTendency_Model(this) + subroutine PreTendency_Model(this) !! PreTendency is a template routine that is used to house any additional calculations !! that you want to execute at the beginning of the tendency calculation routine. !! This default PreTendency simply returns back to the caller without executing any instructions @@ -284,323 +310,323 @@ subroutine PreTendency_Model(this) !! The intention is to provide a method that can be overridden through type-extension, to handle !! any steps that need to be executed before proceeding with the usual tendency calculation methods. !! - implicit none - class(Model),intent(inout) :: this - - return - - endsubroutine PreTendency_Model - - pure function entropy_func_Model(this,s) result(e) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec) :: e - - e = 0.0_prec - - endfunction entropy_func_Model - - pure function riemannflux1d_Model(this,sL,sR,dsdx,nhat) result(flux) - class(Model),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) - real(prec),intent(in) :: nhat - real(prec) :: flux(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - flux(ivar) = 0.0_prec - enddo - - endfunction riemannflux1d_Model - - pure function riemannflux2d_Model(this,sL,sR,dsdx,nhat) result(flux) - class(Model),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 - integer :: ivar - - do ivar = 1,this%nvar - flux(ivar) = 0.0_prec - enddo - - endfunction riemannflux2d_Model - - pure function riemannflux3d_Model(this,sL,sR,dsdx,nhat) result(flux) - class(Model),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:3) - real(prec),intent(in) :: nhat(1:3) - real(prec) :: flux(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - flux(ivar) = 0.0_prec - enddo - - endfunction riemannflux3d_Model - - pure function flux1d_Model(this,s,dsdx) result(flux) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: dsdx(1:this%nvar) - real(prec) :: flux(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - flux(ivar) = 0.0_prec - enddo - - endfunction flux1d_Model - - pure function flux2d_Model(this,s,dsdx) result(flux) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: dsdx(1:this%nvar,1:2) - real(prec) :: flux(1:this%nvar,1:2) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - flux(ivar,1:2) = 0.0_prec - enddo - - endfunction flux2d_Model - - pure function flux3d_Model(this,s,dsdx) result(flux) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: dsdx(1:this%nvar,1:3) - real(prec) :: flux(1:this%nvar,1:3) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - flux(ivar,1:3) = 0.0_prec - enddo - - endfunction flux3d_Model - - pure function source1d_Model(this,s,dsdx) result(source) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: dsdx(1:this%nvar) - real(prec) :: source(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - source(ivar) = 0.0_prec - enddo - - endfunction source1d_Model - - pure function source2d_Model(this,s,dsdx) result(source) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: dsdx(1:this%nvar,1:2) - real(prec) :: source(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - source(ivar) = 0.0_prec - enddo - - endfunction source2d_Model - - pure function source3d_Model(this,s,dsdx) result(source) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: dsdx(1:this%nvar,1:3) - real(prec) :: source(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - source(ivar) = 0.0_prec - enddo - - endfunction source3d_Model - - pure function hbc1d_Generic_Model(this,s,nhat) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - exts(ivar) = 0.0_prec - enddo - - endfunction hbc1d_Generic_Model - - pure function hbc1d_Prescribed_Model(this,x,t) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: x - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - exts(ivar) = 0.0_prec - enddo - - endfunction hbc1d_Prescribed_Model - - pure function hbc2d_Generic_Model(this,s,nhat) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat(1:2) - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - exts(ivar) = 0.0_prec - enddo - - endfunction hbc2d_Generic_Model - - pure function hbc2d_Prescribed_Model(this,x,t) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: x(1:2) - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - exts(ivar) = 0.0_prec - enddo - - endfunction hbc2d_Prescribed_Model - - pure function hbc3d_Generic_Model(this,s,nhat) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat(1:3) - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - exts(ivar) = 0.0_prec - enddo - - endfunction hbc3d_Generic_Model - - pure function hbc3d_Prescribed_Model(this,x,t) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: x(1:3) - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - exts(ivar) = 0.0_prec - enddo - - endfunction hbc3d_Prescribed_Model - - pure function pbc1d_Generic_Model(this,dsdx,nhat) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: dsdx(1:this%nvar) - real(prec),intent(in) :: nhat - real(prec) :: extDsdx(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - extDsdx(ivar) = dsdx(ivar) - enddo - - endfunction pbc1d_Generic_Model - - pure function pbc1d_Prescribed_Model(this,x,t) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: x - real(prec),intent(in) :: t - real(prec) :: extDsdx(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - extDsdx(ivar) = 0.0_prec - enddo - - endfunction pbc1d_Prescribed_Model - - pure function pbc2d_Generic_Model(this,dsdx,nhat) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: dsdx(1:this%nvar,1:2) - real(prec),intent(in) :: nhat(1:2) - real(prec) :: extDsdx(1:this%nvar,1:2) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - extDsdx(ivar,1:2) = dsdx(ivar,1:2) - enddo - - endfunction pbc2d_Generic_Model - - pure function pbc2d_Prescribed_Model(this,x,t) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: x(1:2) - real(prec),intent(in) :: t - real(prec) :: extDsdx(1:this%nvar,1:2) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - extDsdx(ivar,1:2) = 0.0_prec - enddo - - endfunction pbc2d_Prescribed_Model - - pure function pbc3d_Generic_Model(this,dsdx,nhat) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: dsdx(1:this%nvar,1:3) - real(prec),intent(in) :: nhat(1:3) - real(prec) :: extDsdx(1:this%nvar,1:3) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - extDsdx(ivar,1:3) = dsdx(ivar,1:3) - enddo - - endfunction pbc3d_Generic_Model - - pure function pbc3d_Prescribed_Model(this,x,t) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: x(1:3) - real(prec),intent(in) :: t - real(prec) :: extDsdx(1:this%nvar,1:3) - ! Local - integer :: ivar - - do ivar = 1,this%nvar - extDsdx(ivar,1:3) = 0.0_prec - enddo - - endfunction pbc3d_Prescribed_Model - - subroutine SetTimeIntegrator_withChar(this,integrator) + implicit none + class(Model), intent(inout) :: this + + return + + end subroutine PreTendency_Model + + pure function entropy_func_Model(this, s) result(e) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec) :: e + + e = 0.0_prec + + end function entropy_func_Model + + pure function riemannflux1d_Model(this, sL, sR, dsdx, nhat) result(flux) + class(Model), 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) + real(prec), intent(in) :: nhat + real(prec) :: flux(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + flux(ivar) = 0.0_prec + end do + + end function riemannflux1d_Model + + pure function riemannflux2d_Model(this, sL, sR, dsdx, nhat) result(flux) + class(Model), 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 + integer :: ivar + + do ivar = 1, this%nvar + flux(ivar) = 0.0_prec + end do + + end function riemannflux2d_Model + + pure function riemannflux3d_Model(this, sL, sR, dsdx, nhat) result(flux) + class(Model), 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:3) + real(prec), intent(in) :: nhat(1:3) + real(prec) :: flux(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + flux(ivar) = 0.0_prec + end do + + end function riemannflux3d_Model + + pure function flux1d_Model(this, s, dsdx) result(flux) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: dsdx(1:this%nvar) + real(prec) :: flux(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + flux(ivar) = 0.0_prec + end do + + end function flux1d_Model + + pure function flux2d_Model(this, s, dsdx) result(flux) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: dsdx(1:this%nvar, 1:2) + real(prec) :: flux(1:this%nvar, 1:2) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + flux(ivar, 1:2) = 0.0_prec + end do + + end function flux2d_Model + + pure function flux3d_Model(this, s, dsdx) result(flux) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: dsdx(1:this%nvar, 1:3) + real(prec) :: flux(1:this%nvar, 1:3) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + flux(ivar, 1:3) = 0.0_prec + end do + + end function flux3d_Model + + pure function source1d_Model(this, s, dsdx) result(source) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: dsdx(1:this%nvar) + real(prec) :: source(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + source(ivar) = 0.0_prec + end do + + end function source1d_Model + + pure function source2d_Model(this, s, dsdx) result(source) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: dsdx(1:this%nvar, 1:2) + real(prec) :: source(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + source(ivar) = 0.0_prec + end do + + end function source2d_Model + + pure function source3d_Model(this, s, dsdx) result(source) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: dsdx(1:this%nvar, 1:3) + real(prec) :: source(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + source(ivar) = 0.0_prec + end do + + end function source3d_Model + + pure function hbc1d_Generic_Model(this, s, nhat) result(exts) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: nhat + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + exts(ivar) = 0.0_prec + end do + + end function hbc1d_Generic_Model + + pure function hbc1d_Prescribed_Model(this, x, t) result(exts) + class(Model), intent(in) :: this + real(prec), intent(in) :: x + real(prec), intent(in) :: t + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + exts(ivar) = 0.0_prec + end do + + end function hbc1d_Prescribed_Model + + pure function hbc2d_Generic_Model(this, s, nhat) result(exts) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: nhat(1:2) + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + exts(ivar) = 0.0_prec + end do + + end function hbc2d_Generic_Model + + pure function hbc2d_Prescribed_Model(this, x, t) result(exts) + class(Model), intent(in) :: this + real(prec), intent(in) :: x(1:2) + real(prec), intent(in) :: t + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + exts(ivar) = 0.0_prec + end do + + end function hbc2d_Prescribed_Model + + pure function hbc3d_Generic_Model(this, s, nhat) result(exts) + class(Model), intent(in) :: this + real(prec), intent(in) :: s(1:this%nvar) + real(prec), intent(in) :: nhat(1:3) + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + exts(ivar) = 0.0_prec + end do + + end function hbc3d_Generic_Model + + pure function hbc3d_Prescribed_Model(this, x, t) result(exts) + class(Model), intent(in) :: this + real(prec), intent(in) :: x(1:3) + real(prec), intent(in) :: t + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + exts(ivar) = 0.0_prec + end do + + end function hbc3d_Prescribed_Model + + pure function pbc1d_Generic_Model(this, dsdx, nhat) result(extDsdx) + class(Model), intent(in) :: this + real(prec), intent(in) :: dsdx(1:this%nvar) + real(prec), intent(in) :: nhat + real(prec) :: extDsdx(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + extDsdx(ivar) = dsdx(ivar) + end do + + end function pbc1d_Generic_Model + + pure function pbc1d_Prescribed_Model(this, x, t) result(extDsdx) + class(Model), intent(in) :: this + real(prec), intent(in) :: x + real(prec), intent(in) :: t + real(prec) :: extDsdx(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + extDsdx(ivar) = 0.0_prec + end do + + end function pbc1d_Prescribed_Model + + pure function pbc2d_Generic_Model(this, dsdx, nhat) result(extDsdx) + class(Model), intent(in) :: this + real(prec), intent(in) :: dsdx(1:this%nvar, 1:2) + real(prec), intent(in) :: nhat(1:2) + real(prec) :: extDsdx(1:this%nvar, 1:2) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + extDsdx(ivar, 1:2) = dsdx(ivar, 1:2) + end do + + end function pbc2d_Generic_Model + + pure function pbc2d_Prescribed_Model(this, x, t) result(extDsdx) + class(Model), intent(in) :: this + real(prec), intent(in) :: x(1:2) + real(prec), intent(in) :: t + real(prec) :: extDsdx(1:this%nvar, 1:2) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + extDsdx(ivar, 1:2) = 0.0_prec + end do + + end function pbc2d_Prescribed_Model + + pure function pbc3d_Generic_Model(this, dsdx, nhat) result(extDsdx) + class(Model), intent(in) :: this + real(prec), intent(in) :: dsdx(1:this%nvar, 1:3) + real(prec), intent(in) :: nhat(1:3) + real(prec) :: extDsdx(1:this%nvar, 1:3) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + extDsdx(ivar, 1:3) = dsdx(ivar, 1:3) + end do + + end function pbc3d_Generic_Model + + pure function pbc3d_Prescribed_Model(this, x, t) result(extDsdx) + class(Model), intent(in) :: this + real(prec), intent(in) :: x(1:3) + real(prec), intent(in) :: t + real(prec) :: extDsdx(1:this%nvar, 1:3) + ! Local + integer :: ivar + + do ivar = 1, this%nvar + extDsdx(ivar, 1:3) = 0.0_prec + end do + + end function pbc3d_Prescribed_Model + + subroutine SetTimeIntegrator_withChar(this, integrator) !! Sets the time integrator method, using a character input !! !! Valid options for integrator are @@ -612,56 +638,56 @@ subroutine SetTimeIntegrator_withChar(this,integrator) !! !! Note that the character provided is not case-sensitive !! - implicit none - class(Model),intent(inout) :: this - character(*),intent(in) :: integrator - ! Local - character(SELF_INTEGRATOR_LENGTH) :: upperCaseInt + implicit none + class(Model), intent(inout) :: this + character(*), intent(in) :: integrator + ! Local + character(SELF_INTEGRATOR_LENGTH) :: upperCaseInt - upperCaseInt = UpperCase(trim(integrator)) + upperCaseInt = UpperCase(trim(integrator)) - select case(trim(upperCaseInt)) + select case (trim(upperCaseInt)) - case("EULER") - this%timeIntegrator => Euler_timeIntegrator + case ("EULER") + this%timeIntegrator => Euler_timeIntegrator - case("RK2") - this%timeIntegrator => LowStorageRK2_timeIntegrator + case ("RK2") + this%timeIntegrator => LowStorageRK2_timeIntegrator - case("RK3") - this%timeIntegrator => LowStorageRK3_timeIntegrator + case ("RK3") + this%timeIntegrator => LowStorageRK3_timeIntegrator - case("RK4") - this%timeIntegrator => LowStorageRK4_timeIntegrator + case ("RK4") + this%timeIntegrator => LowStorageRK4_timeIntegrator - case DEFAULT - this%timeIntegrator => LowStorageRK3_timeIntegrator + case DEFAULT + this%timeIntegrator => LowStorageRK3_timeIntegrator - endselect + end select - endsubroutine SetTimeIntegrator_withChar + end subroutine SetTimeIntegrator_withChar - subroutine GetSimulationTime(this,t) + subroutine GetSimulationTime(this, t) !! Returns the current simulation time stored in the model % t attribute - implicit none - class(Model),intent(in) :: this - real(prec),intent(out) :: t + implicit none + class(Model), intent(in) :: this + real(prec), intent(out) :: t - t = this%t + t = this%t - endsubroutine GetSimulationTime + end subroutine GetSimulationTime - subroutine SetSimulationTime(this,t) + subroutine SetSimulationTime(this, t) !! Sets the model % t attribute with the provided simulation time - implicit none - class(Model),intent(inout) :: this - real(prec),intent(in) :: t + implicit none + class(Model), intent(inout) :: this + real(prec), intent(in) :: t - this%t = t + this%t = t - endsubroutine SetSimulationTime + end subroutine SetSimulationTime - subroutine CalculateEntropy_Model(this) + subroutine CalculateEntropy_Model(this) !! Base method for calculating entropy of a model !! When this method is not overridden, the entropy !! is simply set to 0.0. When you develop a model @@ -669,45 +695,61 @@ subroutine CalculateEntropy_Model(this) !! children, it is recommended that you define a !! convex mathematical entropy function that is used !! as a measure of the model stability. - implicit none - class(Model),intent(inout) :: this + implicit none + class(Model), intent(inout) :: this - this%entropy = 0.0_prec + this%entropy = 0.0_prec - endsubroutine CalculateEntropy_Model + end subroutine CalculateEntropy_Model - subroutine ReportEntropy_Model(this) + subroutine ReportEntropy_Model(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(Model),intent(in) :: this - ! Local - character(len=20) :: modelTime - character(len=20) :: entropy - character(len=:),allocatable :: str - - ! Copy the time and entropy to a string - write(modelTime,"(ES16.7E3)") this%t - write(entropy,"(ES16.7E3)") this%entropy - - ! Write the output to STDOUT - open(output_unit,ENCODING='utf-8') - write(output_unit,'(A," : ")',ADVANCE='no') __FILE__ - str = 'tᵢ ='//trim(modelTime) - write(output_unit,'(A)',ADVANCE='no') str - str = ' | eᵢ ='//trim(entropy) - write(output_unit,'(A)',ADVANCE='yes') str - - endsubroutine ReportEntropy_Model - - ! ////////////////////////////////////// ! - ! Time Integrators ! - - subroutine ForwardStep_Model(this,tn,dt,ioInterval) + implicit none + class(Model), intent(in) :: this + ! Local + character(len=20) :: modelTime + character(len=20) :: entropy + character(len=:), allocatable :: str + + ! Copy the time and entropy to a string + write (modelTime, "(ES16.7E3)") this%t + write (entropy, "(ES16.7E3)") this%entropy + + ! 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 = ' | eᵢ ='//trim(entropy) + write (output_unit, '(A)', ADVANCE='yes') str + + end subroutine ReportEntropy_Model + + subroutine ReportMetrics_Model(this) + !! Method that can be overridden by users to + !! report their own custom metrics after file io + implicit none + class(Model), intent(inout) :: this + return + end subroutine ReportMetrics_Model + + subroutine ReportUserMetrics_Model(this) + !! Method that can be overridden by users to + !! report their own custom metrics after file io + implicit none + class(Model), intent(inout) :: this + return + end subroutine ReportUserMetrics_Model + + ! ////////////////////////////////////// ! + ! Time Integrators ! + + subroutine ForwardStep_Model(this, tn, dt, ioInterval) !! Forward steps the model using the associated tendency procedure and time integrator !! !! If the final time is provided, the model is forward stepped to that final time, @@ -718,143 +760,149 @@ subroutine ForwardStep_Model(this,tn,dt,ioInterval) !! !! If ioInterval is provided, file IO will be conducted every ioInterval seconds until tn !! is reached - implicit none - class(Model),intent(inout) :: this - real(prec),intent(in) :: tn - real(prec),intent(in) :: dt - real(prec),intent(in) :: ioInterval - ! Local - real(prec) :: targetTime,tNext - integer :: i,nIO - - this%dt = dt - targetTime = tn - - nIO = int((targetTime-this%t)/ioInterval) - do i = 1,nIO - tNext = this%t+ioInterval - call this%timeIntegrator(tNext) - this%t = tNext - call this%WriteModel() - call this%WriteTecplot() - call this%IncrementIOCounter() - call this%CalculateEntropy() - call this%ReportEntropy() - enddo - - endsubroutine ForwardStep_Model - - subroutine Euler_timeIntegrator(this,tn) - implicit none - class(Model),intent(inout) :: this - real(prec),intent(in) :: tn - ! Local - real(prec) :: tRemain - real(prec) :: dtLim - - dtLim = this%dt ! Get the max time step size from the dt attribute - do while(this%t < tn) - - tRemain = tn-this%t - this%dt = min(dtLim,tRemain) - call this%CalculateTendency() - call this%UpdateSolution() - this%t = this%t+this%dt - - enddo - - this%dt = dtLim - - endsubroutine Euler_timeIntegrator - - subroutine LowStorageRK2_timeIntegrator(this,tn) - implicit none - class(Model),intent(inout) :: this - real(prec),intent(in) :: tn - ! Local - integer :: m - real(prec) :: tRemain - real(prec) :: dtLim - real(prec) :: t0 - - dtLim = this%dt ! Get the max time step size from the dt attribute - do while(this%t < tn) - - t0 = this%t - tRemain = tn-this%t - this%dt = min(dtLim,tRemain) - do m = 1,2 - call this%CalculateTendency() - call this%UpdateGRK2(m) - this%t = t0+rk2_b(m)*this%dt - enddo - - this%t = t0+this%dt - - enddo - - this%dt = dtLim - - endsubroutine LowStorageRK2_timeIntegrator - - subroutine LowStorageRK3_timeIntegrator(this,tn) - implicit none - class(Model),intent(inout) :: this - real(prec),intent(in) :: tn - ! Local - integer :: m - real(prec) :: tRemain - real(prec) :: dtLim - real(prec) :: t0 - - dtLim = this%dt ! Get the max time step size from the dt attribute - do while(this%t < tn) - - t0 = this%t - tRemain = tn-this%t - this%dt = min(dtLim,tRemain) - do m = 1,3 - call this%CalculateTendency() - call this%UpdateGRK3(m) - this%t = t0+rk3_b(m)*this%dt - enddo - - this%t = t0+this%dt - - enddo - - this%dt = dtLim - - endsubroutine LowStorageRK3_timeIntegrator - - subroutine LowStorageRK4_timeIntegrator(this,tn) - implicit none - class(Model),intent(inout) :: this - real(prec),intent(in) :: tn - ! Local - integer :: m - real(prec) :: tRemain - real(prec) :: dtLim - real(prec) :: t0 - - dtLim = this%dt ! Get the max time step size from the dt attribute - do while(this%t < tn) + implicit none + class(Model), intent(inout) :: this + real(prec), intent(in) :: tn + real(prec), intent(in) :: dt + real(prec), intent(in) :: ioInterval + ! Local + real(prec) :: targetTime, tNext + integer :: i, nIO + + this%dt = dt + targetTime = tn + + nIO = int((targetTime - this%t)/ioInterval) + do i = 1, nIO + tNext = this%t + ioInterval + call this%timeIntegrator(tNext) + this%t = tNext + + call this%CalculateEntropy() + call this%ReportEntropy() + call this%ReportMetrics() + + call this%WriteModel() + if (this%tecplot_enabled) then + call this%WriteTecplot() + end if + call this%IncrementIOCounter() + + end do + + end subroutine ForwardStep_Model + + subroutine Euler_timeIntegrator(this, tn) + implicit none + class(Model), intent(inout) :: this + real(prec), intent(in) :: tn + ! Local + real(prec) :: tRemain + real(prec) :: dtLim + + dtLim = this%dt ! Get the max time step size from the dt attribute + do while (this%t < tn) + + tRemain = tn - this%t + this%dt = min(dtLim, tRemain) + call this%CalculateTendency() + call this%UpdateSolution() + this%t = this%t + this%dt + + end do + + this%dt = dtLim + + end subroutine Euler_timeIntegrator + + subroutine LowStorageRK2_timeIntegrator(this, tn) + implicit none + class(Model), intent(inout) :: this + real(prec), intent(in) :: tn + ! Local + integer :: m + real(prec) :: tRemain + real(prec) :: dtLim + real(prec) :: t0 + + dtLim = this%dt ! Get the max time step size from the dt attribute + do while (this%t < tn) + + t0 = this%t + tRemain = tn - this%t + this%dt = min(dtLim, tRemain) + do m = 1, 2 + call this%CalculateTendency() + call this%UpdateGRK2(m) + this%t = t0 + rk2_b(m)*this%dt + end do + + this%t = t0 + this%dt + + end do + + this%dt = dtLim + + end subroutine LowStorageRK2_timeIntegrator + + subroutine LowStorageRK3_timeIntegrator(this, tn) + implicit none + class(Model), intent(inout) :: this + real(prec), intent(in) :: tn + ! Local + integer :: m + real(prec) :: tRemain + real(prec) :: dtLim + real(prec) :: t0 + + dtLim = this%dt ! Get the max time step size from the dt attribute + do while (this%t < tn) + + t0 = this%t + tRemain = tn - this%t + this%dt = min(dtLim, tRemain) + do m = 1, 3 + call this%CalculateTendency() + call this%UpdateGRK3(m) + this%t = t0 + rk3_b(m)*this%dt + end do + + this%t = t0 + this%dt + + end do + + this%dt = dtLim + + end subroutine LowStorageRK3_timeIntegrator + + subroutine LowStorageRK4_timeIntegrator(this, tn) + implicit none + class(Model), intent(inout) :: this + real(prec), intent(in) :: tn + ! Local + integer :: m + real(prec) :: tRemain + real(prec) :: dtLim + real(prec) :: t0 + + dtLim = this%dt ! Get the max time step size from the dt attribute + do while (this%t < tn) - t0 = this%t - tRemain = tn-this%t - this%dt = min(dtLim,tRemain) - do m = 1,5 - call this%CalculateTendency() - call this%UpdateGRK4(m) - this%t = t0+rk4_b(m)*this%dt - enddo + t0 = this%t + tRemain = tn - this%t + this%dt = min(dtLim, tRemain) + do m = 1, 5 + call this%CalculateTendency() + call this%UpdateGRK4(m) + this%t = t0 + rk4_b(m)*this%dt + end do - this%t = t0+this%dt + this%t = t0 + this%dt - enddo + end do - this%dt = dtLim + this%dt = dtLim - endsubroutine LowStorageRK4_timeIntegrator + end subroutine LowStorageRK4_timeIntegrator -endmodule SELF_Model +end module SELF_Model From c22e027f2bf3476608e30438e841fecbf24f5d08 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 01:40:34 +0000 Subject: [PATCH 02/11] Fix import bug; add imageio dependency --- pyself/model2d.py | 82 +++++++++++++++++++++++++++++++++-------------- setup.py | 3 +- 2 files changed, 60 insertions(+), 25 deletions(-) diff --git a/pyself/model2d.py b/pyself/model2d.py index 196f8fb85..23a9262c1 100644 --- a/pyself/model2d.py +++ b/pyself/model2d.py @@ -3,7 +3,7 @@ # Other SELF modules -import self.geometry as geometry +import pyself.geometry as geometry class model: def __init__(self): @@ -11,7 +11,6 @@ def __init__(self): 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) diff --git a/setup.py b/setup.py index a17402fea..73efe4628 100644 --- a/setup.py +++ b/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', From 408517219a8be642ad19ea503648124545427c7c Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 01:41:07 +0000 Subject: [PATCH 03/11] Use sphericalsoundwave procedure --- ...ler2d_spherical_soundwave_closeddomain.f90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 b/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 index 098253b2b..a74439c9b 100644 --- a/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 +++ b/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() From a0c9fb9547a93d338e34586a8f7afd76166b358b Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 01:41:30 +0000 Subject: [PATCH 04/11] Disable optimization reports in default flags --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index cd54a6daa..63e3bab50 100644 --- a/CMakeLists.txt +++ b/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" ) From d7a5eb037823d1c5d840edf8c9c6b9bbd4b18239 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 01:48:18 +0000 Subject: [PATCH 05/11] Formatting --- ...ler2d_spherical_soundwave_closeddomain.f90 | 2 +- src/SELF_DGModel2D_t.f90 | 37 +- src/SELF_DGModel3D_t.f90 | 36 +- src/SELF_LinearEuler2D_t.f90 | 50 +- src/SELF_Model.f90 | 1548 ++++++++--------- 5 files changed, 837 insertions(+), 836 deletions(-) diff --git a/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 b/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 index a74439c9b..f2cf58877 100644 --- a/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 +++ b/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 @@ -32,7 +32,7 @@ program LinearEuler_Example 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 :: 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 diff --git a/src/SELF_DGModel2D_t.f90 b/src/SELF_DGModel2D_t.f90 index 394262234..41a44132b 100644 --- a/src/SELF_DGModel2D_t.f90 +++ b/src/SELF_DGModel2D_t.f90 @@ -158,33 +158,34 @@ subroutine ReportMetrics_DGModel2D_t(this) !! open a pull request with modifications to this base !! method. implicit none - class(DGModel2D_t), intent(inout) :: this + class(DGModel2D_t),intent(inout) :: this ! Local character(len=20) :: modelTime - character(len=20) :: minv, maxv - character(len=:), allocatable :: str + 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(modelTime,"(ES16.7E3)") this%t - ! 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 + 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() - end subroutine ReportMetrics_DGModel2D_t + endsubroutine ReportMetrics_DGModel2D_t subroutine SetSolutionFromEqn_DGModel2D_t(this,eqn) implicit none diff --git a/src/SELF_DGModel3D_t.f90 b/src/SELF_DGModel3D_t.f90 index 7d818c1f4..f8c4b13aa 100644 --- a/src/SELF_DGModel3D_t.f90 +++ b/src/SELF_DGModel3D_t.f90 @@ -156,33 +156,33 @@ subroutine ReportMetrics_DGModel3D_t(this) !! open a pull request with modifications to this base !! method. implicit none - class(DGModel3D_t), intent(inout) :: this + class(DGModel3D_t),intent(inout) :: this ! Local character(len=20) :: modelTime - character(len=20) :: minv, maxv - character(len=:), allocatable :: str + 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(modelTime,"(ES16.7E3)") this%t - ! 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 + 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() - end subroutine ReportMetrics_DGModel3D_t + endsubroutine ReportMetrics_DGModel3D_t subroutine SetSolutionFromEqn_DGModel3D_t(this,eqn) implicit none diff --git a/src/SELF_LinearEuler2D_t.f90 b/src/SELF_LinearEuler2D_t.f90 index c5cfdcf1f..3d21835d8 100644 --- a/src/SELF_LinearEuler2D_t.f90 +++ b/src/SELF_LinearEuler2D_t.f90 @@ -189,7 +189,7 @@ 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) + 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 !! @@ -202,36 +202,36 @@ subroutine SphericalSoundWave_LinearEuler2D_t(this, rhoprime, Lr, x0, y0) !! \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 + 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 + 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) + 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) + 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 + 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 + enddo - call this%ReportMetrics() + call this%ReportMetrics() - end subroutine SphericalSoundWave_LinearEuler2D_t + endsubroutine SphericalSoundWave_LinearEuler2D_t endmodule self_LinearEuler2D_t diff --git a/src/SELF_Model.f90 b/src/SELF_Model.f90 index 8fe276bf1..ca8c0bd17 100644 --- a/src/SELF_Model.f90 +++ b/src/SELF_Model.f90 @@ -26,283 +26,283 @@ module SELF_Model - use SELF_SupportRoutines - use SELF_Metadata - use SELF_HDF5 - use HDF5 - use FEQParse + use SELF_SupportRoutines + use SELF_Metadata + use SELF_HDF5 + use HDF5 + use FEQParse - implicit none + implicit none #include "SELF_Macros.h" ! //////////////////////////////////////////////// ! ! Time integration parameters - ! Runge-Kutta 2nd Order (Low Storage) - real(prec), parameter :: rk2_a(1:2) = (/0.0_prec, -0.5_prec/) - real(prec), parameter :: rk2_b(1:2) = (/0.5_prec, 0.5_prec/) - real(prec), parameter :: rk2_g(1:2) = (/0.5_prec, 1.0_prec/) - - ! Williamson's Runge-Kutta 3rd Order (Low Storage) - real(prec), parameter :: rk3_a(1:3) = (/0.0_prec, -5.0_prec/9.0_prec, -153.0_prec/128.0_prec/) - real(prec), parameter :: rk3_b(1:3) = (/0.0_prec, 1.0_prec/3.0_prec, 3.0_prec/4.0_prec/) - real(prec), parameter :: rk3_g(1:3) = (/1.0_prec/3.0_prec, 15.0_prec/16.0_prec, 8.0_prec/15.0_prec/) - - ! Carpenter-Kennedy Runge-Kuttta 4th Order (Low Storage) - real(prec), parameter :: rk4_a(1:5) = (/0.0_prec, & - -1.0_prec, & - -1.0_prec/3.0_prec + & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - -2.0_prec**(1.0_prec/3.0_prec) - & - 2.0_prec**(2.0_prec/3.0_prec) - 2.0_prec, & - -1.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/) - - real(prec), parameter :: rk4_b(1:5) = (/ & - 0.0_prec, & - 2.0_prec/3.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec + & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - 2.0_prec/3.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec + & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - 1.0_prec/3.0_prec - 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec - & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - 1.0_prec/) - - real(prec), parameter :: rk4_g(1:5) = (/ & - 2.0_prec/3.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec + & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - -2.0_prec**(2.0_prec/3.0_prec)/6.0_prec + 1.0_prec/6.0_prec, & - -1.0_prec/3.0_prec - 2.0_prec*2.0_prec**(1.0_prec/3.0_prec)/3.0_prec - & - 2.0_prec**(2.0_prec/3.0_prec)/3.0_prec, & - 1.0_prec/3.0_prec - 2.0_prec**(1.0_prec/3.0_prec)/3.0_prec - & - 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & - 1.0_prec/3.0_prec + 2.0_prec**(1.0_prec/3.0_prec)/6.0_prec + & - 2.0_prec**(2.0_prec/3.0_prec)/12.0_prec/) + ! Runge-Kutta 2nd Order (Low Storage) + real(prec),parameter :: rk2_a(1:2) = (/0.0_prec,-0.5_prec/) + real(prec),parameter :: rk2_b(1:2) = (/0.5_prec,0.5_prec/) + real(prec),parameter :: rk2_g(1:2) = (/0.5_prec,1.0_prec/) + + ! Williamson's Runge-Kutta 3rd Order (Low Storage) + real(prec),parameter :: rk3_a(1:3) = (/0.0_prec,-5.0_prec/9.0_prec,-153.0_prec/128.0_prec/) + real(prec),parameter :: rk3_b(1:3) = (/0.0_prec,1.0_prec/3.0_prec,3.0_prec/4.0_prec/) + real(prec),parameter :: rk3_g(1:3) = (/1.0_prec/3.0_prec,15.0_prec/16.0_prec,8.0_prec/15.0_prec/) + + ! Carpenter-Kennedy Runge-Kuttta 4th Order (Low Storage) + real(prec),parameter :: rk4_a(1:5) = (/0.0_prec, & + -1.0_prec, & + -1.0_prec/3.0_prec+ & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + -2.0_prec**(1.0_prec/3.0_prec)- & + 2.0_prec**(2.0_prec/3.0_prec)-2.0_prec, & + -1.0_prec+2.0_prec**(1.0_prec/3.0_prec)/) + + real(prec),parameter :: rk4_b(1:5) = (/ & + 0.0_prec, & + 2.0_prec/3.0_prec+2.0_prec**(1.0_prec/3.0_prec)/3.0_prec+ & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + 2.0_prec/3.0_prec+2.0_prec**(1.0_prec/3.0_prec)/3.0_prec+ & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + 1.0_prec/3.0_prec-2.0_prec**(1.0_prec/3.0_prec)/3.0_prec- & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + 1.0_prec/) + + real(prec),parameter :: rk4_g(1:5) = (/ & + 2.0_prec/3.0_prec+2.0_prec**(1.0_prec/3.0_prec)/3.0_prec+ & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + -2.0_prec**(2.0_prec/3.0_prec)/6.0_prec+1.0_prec/6.0_prec, & + -1.0_prec/3.0_prec-2.0_prec*2.0_prec**(1.0_prec/3.0_prec)/3.0_prec- & + 2.0_prec**(2.0_prec/3.0_prec)/3.0_prec, & + 1.0_prec/3.0_prec-2.0_prec**(1.0_prec/3.0_prec)/3.0_prec- & + 2.0_prec**(2.0_prec/3.0_prec)/6.0_prec, & + 1.0_prec/3.0_prec+2.0_prec**(1.0_prec/3.0_prec)/6.0_prec+ & + 2.0_prec**(2.0_prec/3.0_prec)/12.0_prec/) ! - integer, parameter :: SELF_EULER = 100 - integer, parameter :: SELF_RK2 = 200 - integer, parameter :: SELF_RK3 = 300 - integer, parameter :: SELF_RK4 = 400 - ! integer,parameter :: SELF_AB2 = 201 - ! integer,parameter :: SELF_AB3 = 301 - ! integer,parameter :: SELF_AB4 = 401 + integer,parameter :: SELF_EULER = 100 + integer,parameter :: SELF_RK2 = 200 + integer,parameter :: SELF_RK3 = 300 + integer,parameter :: SELF_RK4 = 400 + ! integer,parameter :: SELF_AB2 = 201 + ! integer,parameter :: SELF_AB3 = 301 + ! integer,parameter :: SELF_AB4 = 401 - integer, parameter :: SELF_INTEGRATOR_LENGTH = 10 ! max length of integrator methods when specified as char - integer, parameter :: SELF_EQUATION_LENGTH = 500 + integer,parameter :: SELF_INTEGRATOR_LENGTH = 10 ! max length of integrator methods when specified as char + integer,parameter :: SELF_EQUATION_LENGTH = 500 ! //////////////////////////////////////////////// ! ! Model Formulations ! - integer, parameter :: SELF_FORMULATION_LENGTH = 30 ! max length of integrator methods when specified as char - - type, abstract :: Model - - ! Time integration attributes - procedure(SELF_timeIntegrator), pointer :: timeIntegrator => Euler_timeIntegrator - real(prec) :: dt - real(prec) :: t - integer :: ioIterate = 0 - logical :: gradient_enabled = .false. - logical :: prescribed_bcs_enabled = .true. - logical :: tecplot_enabled = .true. - integer :: nvar - ! Standard Diagnostics - real(prec) :: entropy ! Mathematical entropy function for the model - - contains - - procedure :: IncrementIOCounter - - procedure :: PrintType => PrintType_Model - - procedure :: SetNumberOfVariables => SetNumberOfVariables_Model - - procedure :: AdditionalInit => AdditionalInit_Model - procedure :: AdditionalFree => AdditionalFree_Model - procedure :: AdditionalOutput => AdditionalOutput_Model - - procedure :: ForwardStep => ForwardStep_Model - - procedure :: Euler_timeIntegrator - - ! Runge-Kutta methods - procedure :: LowStorageRK2_timeIntegrator - procedure(UpdateGRK), deferred :: UpdateGRK2 - - procedure :: LowStorageRK3_timeIntegrator - procedure(UpdateGRK), deferred :: UpdateGRK3 - - procedure :: LowStorageRK4_timeIntegrator - procedure(UpdateGRK), deferred :: UpdateGRK4 - - procedure :: PreTendency => PreTendency_Model - procedure :: entropy_func => entropy_func_Model - - procedure :: flux1D => flux1d_Model - procedure :: flux2D => flux2d_Model - procedure :: flux3D => flux3d_Model - - procedure :: riemannflux1d => riemannflux1d_Model - procedure :: riemannflux2d => riemannflux2d_Model - procedure :: riemannflux3d => riemannflux3d_Model - - procedure :: source1d => source1d_Model - procedure :: source2d => source2d_Model - procedure :: source3d => source3d_Model - - ! Boundary condition functions (hyperbolic) - procedure :: hbc1d_Prescribed => hbc1d_Prescribed_Model - procedure :: hbc1d_Radiation => hbc1d_Generic_Model - procedure :: hbc1d_NoNormalFlow => hbc1d_Generic_Model - procedure :: hbc2d_Prescribed => hbc2d_Prescribed_Model - procedure :: hbc2d_Radiation => hbc2d_Generic_Model - procedure :: hbc2d_NoNormalFlow => hbc2d_Generic_Model - procedure :: hbc3d_Prescribed => hbc3d_Prescribed_Model - procedure :: hbc3d_Radiation => hbc3d_Generic_Model - procedure :: hbc3d_NoNormalFlow => hbc3d_Generic_Model - - ! Boundary condition functions (parabolic) - procedure :: pbc1d_Prescribed => pbc1d_Prescribed_Model - procedure :: pbc1d_Radiation => pbc1d_Generic_Model - procedure :: pbc1d_NoNormalFlow => pbc1d_Generic_Model - procedure :: pbc2d_Prescribed => pbc2d_Prescribed_Model - procedure :: pbc2d_Radiation => pbc2d_Generic_Model - procedure :: pbc2d_NoNormalFlow => pbc2d_Generic_Model - procedure :: pbc3d_Prescribed => pbc3d_Prescribed_Model - procedure :: pbc3d_Radiation => pbc3d_Generic_Model - procedure :: pbc3d_NoNormalFlow => pbc3d_Generic_Model - - procedure :: ReportEntropy => ReportEntropy_Model - procedure :: ReportMetrics => ReportMetrics_Model - procedure :: ReportUserMetrics => ReportUserMetrics_Model - procedure :: CalculateEntropy => CalculateEntropy_Model - - procedure(UpdateSolution), deferred :: UpdateSolution - procedure(CalculateTendency), deferred :: CalculateTendency - procedure(ReadModel), deferred :: ReadModel - procedure(WriteModel), deferred :: WriteModel - procedure(WriteTecplot), deferred :: WriteTecplot - - generic :: SetTimeIntegrator => SetTimeIntegrator_withChar - procedure, private :: SetTimeIntegrator_withChar - - procedure :: SetSimulationTime - procedure :: GetSimulationTime - - end type Model - - interface - subroutine SELF_timeIntegrator(this, tn) - use SELF_Constants, only: prec - import Model - implicit none - class(Model), intent(inout) :: this - real(prec), intent(in) :: tn - end subroutine SELF_timeIntegrator - end interface - - interface - subroutine UpdateGRK(this, m) - import Model - implicit none - class(Model), intent(inout) :: this - integer, intent(in) :: m - end subroutine UpdateGRK - end interface - - interface - subroutine UpdateSolution(this, dt) - use SELF_Constants, only: prec - import Model - implicit none - class(Model), intent(inout) :: this - real(prec), optional, intent(in) :: dt - end subroutine UpdateSolution - end interface - - interface - subroutine CalculateTendency(this) - import Model - implicit none - class(Model), intent(inout) :: this - end subroutine CalculateTendency - end interface - - interface - subroutine WriteModel(this, filename) - import Model - implicit none - class(Model), intent(inout) :: this - character(*), intent(in), optional :: filename - end subroutine WriteModel - end interface - - interface - subroutine ReadModel(this, filename) - import Model - implicit none - class(Model), intent(inout) :: this - character(*), intent(in) :: filename - end subroutine ReadModel - end interface - - interface - subroutine WriteTecplot(this, filename) - import Model - implicit none - class(Model), intent(inout) :: this - character(*), intent(in), optional :: filename - end subroutine WriteTecplot - end interface + integer,parameter :: SELF_FORMULATION_LENGTH = 30 ! max length of integrator methods when specified as char + + type,abstract :: Model + + ! Time integration attributes + procedure(SELF_timeIntegrator),pointer :: timeIntegrator => Euler_timeIntegrator + real(prec) :: dt + real(prec) :: t + integer :: ioIterate = 0 + logical :: gradient_enabled = .false. + logical :: prescribed_bcs_enabled = .true. + logical :: tecplot_enabled = .true. + integer :: nvar + ! Standard Diagnostics + real(prec) :: entropy ! Mathematical entropy function for the model + + contains + + procedure :: IncrementIOCounter + + procedure :: PrintType => PrintType_Model + + procedure :: SetNumberOfVariables => SetNumberOfVariables_Model + + procedure :: AdditionalInit => AdditionalInit_Model + procedure :: AdditionalFree => AdditionalFree_Model + procedure :: AdditionalOutput => AdditionalOutput_Model + + procedure :: ForwardStep => ForwardStep_Model + + procedure :: Euler_timeIntegrator + + ! Runge-Kutta methods + procedure :: LowStorageRK2_timeIntegrator + procedure(UpdateGRK),deferred :: UpdateGRK2 + + procedure :: LowStorageRK3_timeIntegrator + procedure(UpdateGRK),deferred :: UpdateGRK3 + + procedure :: LowStorageRK4_timeIntegrator + procedure(UpdateGRK),deferred :: UpdateGRK4 + + procedure :: PreTendency => PreTendency_Model + procedure :: entropy_func => entropy_func_Model + + procedure :: flux1D => flux1d_Model + procedure :: flux2D => flux2d_Model + procedure :: flux3D => flux3d_Model + + procedure :: riemannflux1d => riemannflux1d_Model + procedure :: riemannflux2d => riemannflux2d_Model + procedure :: riemannflux3d => riemannflux3d_Model + + procedure :: source1d => source1d_Model + procedure :: source2d => source2d_Model + procedure :: source3d => source3d_Model + + ! Boundary condition functions (hyperbolic) + procedure :: hbc1d_Prescribed => hbc1d_Prescribed_Model + procedure :: hbc1d_Radiation => hbc1d_Generic_Model + procedure :: hbc1d_NoNormalFlow => hbc1d_Generic_Model + procedure :: hbc2d_Prescribed => hbc2d_Prescribed_Model + procedure :: hbc2d_Radiation => hbc2d_Generic_Model + procedure :: hbc2d_NoNormalFlow => hbc2d_Generic_Model + procedure :: hbc3d_Prescribed => hbc3d_Prescribed_Model + procedure :: hbc3d_Radiation => hbc3d_Generic_Model + procedure :: hbc3d_NoNormalFlow => hbc3d_Generic_Model + + ! Boundary condition functions (parabolic) + procedure :: pbc1d_Prescribed => pbc1d_Prescribed_Model + procedure :: pbc1d_Radiation => pbc1d_Generic_Model + procedure :: pbc1d_NoNormalFlow => pbc1d_Generic_Model + procedure :: pbc2d_Prescribed => pbc2d_Prescribed_Model + procedure :: pbc2d_Radiation => pbc2d_Generic_Model + procedure :: pbc2d_NoNormalFlow => pbc2d_Generic_Model + procedure :: pbc3d_Prescribed => pbc3d_Prescribed_Model + procedure :: pbc3d_Radiation => pbc3d_Generic_Model + procedure :: pbc3d_NoNormalFlow => pbc3d_Generic_Model + + procedure :: ReportEntropy => ReportEntropy_Model + procedure :: ReportMetrics => ReportMetrics_Model + procedure :: ReportUserMetrics => ReportUserMetrics_Model + procedure :: CalculateEntropy => CalculateEntropy_Model + + procedure(UpdateSolution),deferred :: UpdateSolution + procedure(CalculateTendency),deferred :: CalculateTendency + procedure(ReadModel),deferred :: ReadModel + procedure(WriteModel),deferred :: WriteModel + procedure(WriteTecplot),deferred :: WriteTecplot + + generic :: SetTimeIntegrator => SetTimeIntegrator_withChar + procedure,private :: SetTimeIntegrator_withChar + + procedure :: SetSimulationTime + procedure :: GetSimulationTime + + endtype Model + + interface + subroutine SELF_timeIntegrator(this,tn) + use SELF_Constants,only:prec + import Model + implicit none + class(Model),intent(inout) :: this + real(prec),intent(in) :: tn + endsubroutine SELF_timeIntegrator + endinterface + + interface + subroutine UpdateGRK(this,m) + import Model + implicit none + class(Model),intent(inout) :: this + integer,intent(in) :: m + endsubroutine UpdateGRK + endinterface + + interface + subroutine UpdateSolution(this,dt) + use SELF_Constants,only:prec + import Model + implicit none + class(Model),intent(inout) :: this + real(prec),optional,intent(in) :: dt + endsubroutine UpdateSolution + endinterface + + interface + subroutine CalculateTendency(this) + import Model + implicit none + class(Model),intent(inout) :: this + endsubroutine CalculateTendency + endinterface + + interface + subroutine WriteModel(this,filename) + import Model + implicit none + class(Model),intent(inout) :: this + character(*),intent(in),optional :: filename + endsubroutine WriteModel + endinterface + + interface + subroutine ReadModel(this,filename) + import Model + implicit none + class(Model),intent(inout) :: this + character(*),intent(in) :: filename + endsubroutine ReadModel + endinterface + + interface + subroutine WriteTecplot(this,filename) + import Model + implicit none + class(Model),intent(inout) :: this + character(*),intent(in),optional :: filename + endsubroutine WriteTecplot + endinterface contains - subroutine IncrementIOCounter(this) - implicit none - class(Model), intent(inout) :: this + subroutine IncrementIOCounter(this) + implicit none + class(Model),intent(inout) :: this - ! Increment the ioIterate - this%ioIterate = this%ioIterate + 1 + ! Increment the ioIterate + this%ioIterate = this%ioIterate+1 - end subroutine IncrementIOCounter + endsubroutine IncrementIOCounter - subroutine SetNumberOfVariables_Model(this) - implicit none - class(Model), intent(inout) :: this + subroutine SetNumberOfVariables_Model(this) + implicit none + class(Model),intent(inout) :: this - this%nvar = 1 + this%nvar = 1 - end subroutine SetNumberOfVariables_Model + endsubroutine SetNumberOfVariables_Model - subroutine AdditionalInit_Model(this) - implicit none - class(Model), intent(inout) :: this - return - end subroutine AdditionalInit_Model + subroutine AdditionalInit_Model(this) + implicit none + class(Model),intent(inout) :: this + return + endsubroutine AdditionalInit_Model - subroutine AdditionalFree_Model(this) - implicit none - class(Model), intent(inout) :: this - return - end subroutine AdditionalFree_Model + subroutine AdditionalFree_Model(this) + implicit none + class(Model),intent(inout) :: this + return + endsubroutine AdditionalFree_Model - subroutine AdditionalOutput_Model(this, fileid) - implicit none - class(Model), intent(inout) :: this - integer(HID_T), intent(in) :: fileid - return - end subroutine AdditionalOutput_Model + subroutine AdditionalOutput_Model(this,fileid) + implicit none + class(Model),intent(inout) :: this + integer(HID_T),intent(in) :: fileid + return + endsubroutine AdditionalOutput_Model - subroutine PrintType_Model(this) - implicit none - class(Model), intent(in) :: this + subroutine PrintType_Model(this) + implicit none + class(Model),intent(in) :: this - print *, __FILE__//" : Model : No model type" + print*,__FILE__//" : Model : No model type" - end subroutine PrintType_Model + endsubroutine PrintType_Model - subroutine PreTendency_Model(this) + subroutine PreTendency_Model(this) !! PreTendency is a template routine that is used to house any additional calculations !! that you want to execute at the beginning of the tendency calculation routine. !! This default PreTendency simply returns back to the caller without executing any instructions @@ -310,323 +310,323 @@ subroutine PreTendency_Model(this) !! The intention is to provide a method that can be overridden through type-extension, to handle !! any steps that need to be executed before proceeding with the usual tendency calculation methods. !! - implicit none - class(Model), intent(inout) :: this - - return - - end subroutine PreTendency_Model - - pure function entropy_func_Model(this, s) result(e) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec) :: e - - e = 0.0_prec - - end function entropy_func_Model - - pure function riemannflux1d_Model(this, sL, sR, dsdx, nhat) result(flux) - class(Model), 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) - real(prec), intent(in) :: nhat - real(prec) :: flux(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - flux(ivar) = 0.0_prec - end do - - end function riemannflux1d_Model - - pure function riemannflux2d_Model(this, sL, sR, dsdx, nhat) result(flux) - class(Model), 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 - integer :: ivar - - do ivar = 1, this%nvar - flux(ivar) = 0.0_prec - end do - - end function riemannflux2d_Model - - pure function riemannflux3d_Model(this, sL, sR, dsdx, nhat) result(flux) - class(Model), 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:3) - real(prec), intent(in) :: nhat(1:3) - real(prec) :: flux(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - flux(ivar) = 0.0_prec - end do - - end function riemannflux3d_Model - - pure function flux1d_Model(this, s, dsdx) result(flux) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: dsdx(1:this%nvar) - real(prec) :: flux(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - flux(ivar) = 0.0_prec - end do - - end function flux1d_Model - - pure function flux2d_Model(this, s, dsdx) result(flux) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: dsdx(1:this%nvar, 1:2) - real(prec) :: flux(1:this%nvar, 1:2) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - flux(ivar, 1:2) = 0.0_prec - end do - - end function flux2d_Model - - pure function flux3d_Model(this, s, dsdx) result(flux) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: dsdx(1:this%nvar, 1:3) - real(prec) :: flux(1:this%nvar, 1:3) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - flux(ivar, 1:3) = 0.0_prec - end do - - end function flux3d_Model - - pure function source1d_Model(this, s, dsdx) result(source) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: dsdx(1:this%nvar) - real(prec) :: source(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - source(ivar) = 0.0_prec - end do - - end function source1d_Model - - pure function source2d_Model(this, s, dsdx) result(source) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: dsdx(1:this%nvar, 1:2) - real(prec) :: source(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - source(ivar) = 0.0_prec - end do - - end function source2d_Model - - pure function source3d_Model(this, s, dsdx) result(source) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: dsdx(1:this%nvar, 1:3) - real(prec) :: source(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - source(ivar) = 0.0_prec - end do - - end function source3d_Model - - pure function hbc1d_Generic_Model(this, s, nhat) result(exts) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: nhat - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - exts(ivar) = 0.0_prec - end do - - end function hbc1d_Generic_Model - - pure function hbc1d_Prescribed_Model(this, x, t) result(exts) - class(Model), intent(in) :: this - real(prec), intent(in) :: x - real(prec), intent(in) :: t - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - exts(ivar) = 0.0_prec - end do - - end function hbc1d_Prescribed_Model - - pure function hbc2d_Generic_Model(this, s, nhat) result(exts) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: nhat(1:2) - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - exts(ivar) = 0.0_prec - end do - - end function hbc2d_Generic_Model - - pure function hbc2d_Prescribed_Model(this, x, t) result(exts) - class(Model), intent(in) :: this - real(prec), intent(in) :: x(1:2) - real(prec), intent(in) :: t - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - exts(ivar) = 0.0_prec - end do - - end function hbc2d_Prescribed_Model - - pure function hbc3d_Generic_Model(this, s, nhat) result(exts) - class(Model), intent(in) :: this - real(prec), intent(in) :: s(1:this%nvar) - real(prec), intent(in) :: nhat(1:3) - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - exts(ivar) = 0.0_prec - end do - - end function hbc3d_Generic_Model - - pure function hbc3d_Prescribed_Model(this, x, t) result(exts) - class(Model), intent(in) :: this - real(prec), intent(in) :: x(1:3) - real(prec), intent(in) :: t - real(prec) :: exts(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - exts(ivar) = 0.0_prec - end do - - end function hbc3d_Prescribed_Model - - pure function pbc1d_Generic_Model(this, dsdx, nhat) result(extDsdx) - class(Model), intent(in) :: this - real(prec), intent(in) :: dsdx(1:this%nvar) - real(prec), intent(in) :: nhat - real(prec) :: extDsdx(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - extDsdx(ivar) = dsdx(ivar) - end do - - end function pbc1d_Generic_Model - - pure function pbc1d_Prescribed_Model(this, x, t) result(extDsdx) - class(Model), intent(in) :: this - real(prec), intent(in) :: x - real(prec), intent(in) :: t - real(prec) :: extDsdx(1:this%nvar) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - extDsdx(ivar) = 0.0_prec - end do - - end function pbc1d_Prescribed_Model - - pure function pbc2d_Generic_Model(this, dsdx, nhat) result(extDsdx) - class(Model), intent(in) :: this - real(prec), intent(in) :: dsdx(1:this%nvar, 1:2) - real(prec), intent(in) :: nhat(1:2) - real(prec) :: extDsdx(1:this%nvar, 1:2) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - extDsdx(ivar, 1:2) = dsdx(ivar, 1:2) - end do - - end function pbc2d_Generic_Model - - pure function pbc2d_Prescribed_Model(this, x, t) result(extDsdx) - class(Model), intent(in) :: this - real(prec), intent(in) :: x(1:2) - real(prec), intent(in) :: t - real(prec) :: extDsdx(1:this%nvar, 1:2) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - extDsdx(ivar, 1:2) = 0.0_prec - end do - - end function pbc2d_Prescribed_Model - - pure function pbc3d_Generic_Model(this, dsdx, nhat) result(extDsdx) - class(Model), intent(in) :: this - real(prec), intent(in) :: dsdx(1:this%nvar, 1:3) - real(prec), intent(in) :: nhat(1:3) - real(prec) :: extDsdx(1:this%nvar, 1:3) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - extDsdx(ivar, 1:3) = dsdx(ivar, 1:3) - end do - - end function pbc3d_Generic_Model - - pure function pbc3d_Prescribed_Model(this, x, t) result(extDsdx) - class(Model), intent(in) :: this - real(prec), intent(in) :: x(1:3) - real(prec), intent(in) :: t - real(prec) :: extDsdx(1:this%nvar, 1:3) - ! Local - integer :: ivar - - do ivar = 1, this%nvar - extDsdx(ivar, 1:3) = 0.0_prec - end do - - end function pbc3d_Prescribed_Model - - subroutine SetTimeIntegrator_withChar(this, integrator) + implicit none + class(Model),intent(inout) :: this + + return + + endsubroutine PreTendency_Model + + pure function entropy_func_Model(this,s) result(e) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec) :: e + + e = 0.0_prec + + endfunction entropy_func_Model + + pure function riemannflux1d_Model(this,sL,sR,dsdx,nhat) result(flux) + class(Model),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) + real(prec),intent(in) :: nhat + real(prec) :: flux(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + flux(ivar) = 0.0_prec + enddo + + endfunction riemannflux1d_Model + + pure function riemannflux2d_Model(this,sL,sR,dsdx,nhat) result(flux) + class(Model),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 + integer :: ivar + + do ivar = 1,this%nvar + flux(ivar) = 0.0_prec + enddo + + endfunction riemannflux2d_Model + + pure function riemannflux3d_Model(this,sL,sR,dsdx,nhat) result(flux) + class(Model),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:3) + real(prec),intent(in) :: nhat(1:3) + real(prec) :: flux(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + flux(ivar) = 0.0_prec + enddo + + endfunction riemannflux3d_Model + + pure function flux1d_Model(this,s,dsdx) result(flux) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: dsdx(1:this%nvar) + real(prec) :: flux(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + flux(ivar) = 0.0_prec + enddo + + endfunction flux1d_Model + + pure function flux2d_Model(this,s,dsdx) result(flux) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: dsdx(1:this%nvar,1:2) + real(prec) :: flux(1:this%nvar,1:2) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + flux(ivar,1:2) = 0.0_prec + enddo + + endfunction flux2d_Model + + pure function flux3d_Model(this,s,dsdx) result(flux) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: dsdx(1:this%nvar,1:3) + real(prec) :: flux(1:this%nvar,1:3) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + flux(ivar,1:3) = 0.0_prec + enddo + + endfunction flux3d_Model + + pure function source1d_Model(this,s,dsdx) result(source) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: dsdx(1:this%nvar) + real(prec) :: source(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + source(ivar) = 0.0_prec + enddo + + endfunction source1d_Model + + pure function source2d_Model(this,s,dsdx) result(source) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: dsdx(1:this%nvar,1:2) + real(prec) :: source(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + source(ivar) = 0.0_prec + enddo + + endfunction source2d_Model + + pure function source3d_Model(this,s,dsdx) result(source) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: dsdx(1:this%nvar,1:3) + real(prec) :: source(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + source(ivar) = 0.0_prec + enddo + + endfunction source3d_Model + + pure function hbc1d_Generic_Model(this,s,nhat) result(exts) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: nhat + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + exts(ivar) = 0.0_prec + enddo + + endfunction hbc1d_Generic_Model + + pure function hbc1d_Prescribed_Model(this,x,t) result(exts) + class(Model),intent(in) :: this + real(prec),intent(in) :: x + real(prec),intent(in) :: t + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + exts(ivar) = 0.0_prec + enddo + + endfunction hbc1d_Prescribed_Model + + pure function hbc2d_Generic_Model(this,s,nhat) result(exts) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: nhat(1:2) + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + exts(ivar) = 0.0_prec + enddo + + endfunction hbc2d_Generic_Model + + pure function hbc2d_Prescribed_Model(this,x,t) result(exts) + class(Model),intent(in) :: this + real(prec),intent(in) :: x(1:2) + real(prec),intent(in) :: t + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + exts(ivar) = 0.0_prec + enddo + + endfunction hbc2d_Prescribed_Model + + pure function hbc3d_Generic_Model(this,s,nhat) result(exts) + class(Model),intent(in) :: this + real(prec),intent(in) :: s(1:this%nvar) + real(prec),intent(in) :: nhat(1:3) + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + exts(ivar) = 0.0_prec + enddo + + endfunction hbc3d_Generic_Model + + pure function hbc3d_Prescribed_Model(this,x,t) result(exts) + class(Model),intent(in) :: this + real(prec),intent(in) :: x(1:3) + real(prec),intent(in) :: t + real(prec) :: exts(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + exts(ivar) = 0.0_prec + enddo + + endfunction hbc3d_Prescribed_Model + + pure function pbc1d_Generic_Model(this,dsdx,nhat) result(extDsdx) + class(Model),intent(in) :: this + real(prec),intent(in) :: dsdx(1:this%nvar) + real(prec),intent(in) :: nhat + real(prec) :: extDsdx(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + extDsdx(ivar) = dsdx(ivar) + enddo + + endfunction pbc1d_Generic_Model + + pure function pbc1d_Prescribed_Model(this,x,t) result(extDsdx) + class(Model),intent(in) :: this + real(prec),intent(in) :: x + real(prec),intent(in) :: t + real(prec) :: extDsdx(1:this%nvar) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + extDsdx(ivar) = 0.0_prec + enddo + + endfunction pbc1d_Prescribed_Model + + pure function pbc2d_Generic_Model(this,dsdx,nhat) result(extDsdx) + class(Model),intent(in) :: this + real(prec),intent(in) :: dsdx(1:this%nvar,1:2) + real(prec),intent(in) :: nhat(1:2) + real(prec) :: extDsdx(1:this%nvar,1:2) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + extDsdx(ivar,1:2) = dsdx(ivar,1:2) + enddo + + endfunction pbc2d_Generic_Model + + pure function pbc2d_Prescribed_Model(this,x,t) result(extDsdx) + class(Model),intent(in) :: this + real(prec),intent(in) :: x(1:2) + real(prec),intent(in) :: t + real(prec) :: extDsdx(1:this%nvar,1:2) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + extDsdx(ivar,1:2) = 0.0_prec + enddo + + endfunction pbc2d_Prescribed_Model + + pure function pbc3d_Generic_Model(this,dsdx,nhat) result(extDsdx) + class(Model),intent(in) :: this + real(prec),intent(in) :: dsdx(1:this%nvar,1:3) + real(prec),intent(in) :: nhat(1:3) + real(prec) :: extDsdx(1:this%nvar,1:3) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + extDsdx(ivar,1:3) = dsdx(ivar,1:3) + enddo + + endfunction pbc3d_Generic_Model + + pure function pbc3d_Prescribed_Model(this,x,t) result(extDsdx) + class(Model),intent(in) :: this + real(prec),intent(in) :: x(1:3) + real(prec),intent(in) :: t + real(prec) :: extDsdx(1:this%nvar,1:3) + ! Local + integer :: ivar + + do ivar = 1,this%nvar + extDsdx(ivar,1:3) = 0.0_prec + enddo + + endfunction pbc3d_Prescribed_Model + + subroutine SetTimeIntegrator_withChar(this,integrator) !! Sets the time integrator method, using a character input !! !! Valid options for integrator are @@ -638,56 +638,56 @@ subroutine SetTimeIntegrator_withChar(this, integrator) !! !! Note that the character provided is not case-sensitive !! - implicit none - class(Model), intent(inout) :: this - character(*), intent(in) :: integrator - ! Local - character(SELF_INTEGRATOR_LENGTH) :: upperCaseInt + implicit none + class(Model),intent(inout) :: this + character(*),intent(in) :: integrator + ! Local + character(SELF_INTEGRATOR_LENGTH) :: upperCaseInt - upperCaseInt = UpperCase(trim(integrator)) + upperCaseInt = UpperCase(trim(integrator)) - select case (trim(upperCaseInt)) + select case(trim(upperCaseInt)) - case ("EULER") - this%timeIntegrator => Euler_timeIntegrator + case("EULER") + this%timeIntegrator => Euler_timeIntegrator - case ("RK2") - this%timeIntegrator => LowStorageRK2_timeIntegrator + case("RK2") + this%timeIntegrator => LowStorageRK2_timeIntegrator - case ("RK3") - this%timeIntegrator => LowStorageRK3_timeIntegrator + case("RK3") + this%timeIntegrator => LowStorageRK3_timeIntegrator - case ("RK4") - this%timeIntegrator => LowStorageRK4_timeIntegrator + case("RK4") + this%timeIntegrator => LowStorageRK4_timeIntegrator - case DEFAULT - this%timeIntegrator => LowStorageRK3_timeIntegrator + case DEFAULT + this%timeIntegrator => LowStorageRK3_timeIntegrator - end select + endselect - end subroutine SetTimeIntegrator_withChar + endsubroutine SetTimeIntegrator_withChar - subroutine GetSimulationTime(this, t) + subroutine GetSimulationTime(this,t) !! Returns the current simulation time stored in the model % t attribute - implicit none - class(Model), intent(in) :: this - real(prec), intent(out) :: t + implicit none + class(Model),intent(in) :: this + real(prec),intent(out) :: t - t = this%t + t = this%t - end subroutine GetSimulationTime + endsubroutine GetSimulationTime - subroutine SetSimulationTime(this, t) + subroutine SetSimulationTime(this,t) !! Sets the model % t attribute with the provided simulation time - implicit none - class(Model), intent(inout) :: this - real(prec), intent(in) :: t + implicit none + class(Model),intent(inout) :: this + real(prec),intent(in) :: t - this%t = t + this%t = t - end subroutine SetSimulationTime + endsubroutine SetSimulationTime - subroutine CalculateEntropy_Model(this) + subroutine CalculateEntropy_Model(this) !! Base method for calculating entropy of a model !! When this method is not overridden, the entropy !! is simply set to 0.0. When you develop a model @@ -695,61 +695,61 @@ subroutine CalculateEntropy_Model(this) !! children, it is recommended that you define a !! convex mathematical entropy function that is used !! as a measure of the model stability. - implicit none - class(Model), intent(inout) :: this + implicit none + class(Model),intent(inout) :: this - this%entropy = 0.0_prec + this%entropy = 0.0_prec - end subroutine CalculateEntropy_Model + endsubroutine CalculateEntropy_Model - subroutine ReportEntropy_Model(this) + subroutine ReportEntropy_Model(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(Model), intent(in) :: this - ! Local - character(len=20) :: modelTime - character(len=20) :: entropy - character(len=:), allocatable :: str - - ! Copy the time and entropy to a string - write (modelTime, "(ES16.7E3)") this%t - write (entropy, "(ES16.7E3)") this%entropy - - ! 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 = ' | eᵢ ='//trim(entropy) - write (output_unit, '(A)', ADVANCE='yes') str - - end subroutine ReportEntropy_Model - - subroutine ReportMetrics_Model(this) + implicit none + class(Model),intent(in) :: this + ! Local + character(len=20) :: modelTime + character(len=20) :: entropy + character(len=:),allocatable :: str + + ! Copy the time and entropy to a string + write(modelTime,"(ES16.7E3)") this%t + write(entropy,"(ES16.7E3)") this%entropy + + ! 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 = ' | eᵢ ='//trim(entropy) + write(output_unit,'(A)',ADVANCE='yes') str + + endsubroutine ReportEntropy_Model + + subroutine ReportMetrics_Model(this) !! Method that can be overridden by users to !! report their own custom metrics after file io - implicit none - class(Model), intent(inout) :: this - return - end subroutine ReportMetrics_Model + implicit none + class(Model),intent(inout) :: this + return + endsubroutine ReportMetrics_Model - subroutine ReportUserMetrics_Model(this) + subroutine ReportUserMetrics_Model(this) !! Method that can be overridden by users to !! report their own custom metrics after file io - implicit none - class(Model), intent(inout) :: this - return - end subroutine ReportUserMetrics_Model + implicit none + class(Model),intent(inout) :: this + return + endsubroutine ReportUserMetrics_Model - ! ////////////////////////////////////// ! - ! Time Integrators ! + ! ////////////////////////////////////// ! + ! Time Integrators ! - subroutine ForwardStep_Model(this, tn, dt, ioInterval) + subroutine ForwardStep_Model(this,tn,dt,ioInterval) !! Forward steps the model using the associated tendency procedure and time integrator !! !! If the final time is provided, the model is forward stepped to that final time, @@ -760,149 +760,149 @@ subroutine ForwardStep_Model(this, tn, dt, ioInterval) !! !! If ioInterval is provided, file IO will be conducted every ioInterval seconds until tn !! is reached - implicit none - class(Model), intent(inout) :: this - real(prec), intent(in) :: tn - real(prec), intent(in) :: dt - real(prec), intent(in) :: ioInterval - ! Local - real(prec) :: targetTime, tNext - integer :: i, nIO - - this%dt = dt - targetTime = tn - - nIO = int((targetTime - this%t)/ioInterval) - do i = 1, nIO - tNext = this%t + ioInterval - call this%timeIntegrator(tNext) - this%t = tNext - - call this%CalculateEntropy() - call this%ReportEntropy() - call this%ReportMetrics() - - call this%WriteModel() - if (this%tecplot_enabled) then - call this%WriteTecplot() - end if - call this%IncrementIOCounter() - - end do - - end subroutine ForwardStep_Model - - subroutine Euler_timeIntegrator(this, tn) - implicit none - class(Model), intent(inout) :: this - real(prec), intent(in) :: tn - ! Local - real(prec) :: tRemain - real(prec) :: dtLim + implicit none + class(Model),intent(inout) :: this + real(prec),intent(in) :: tn + real(prec),intent(in) :: dt + real(prec),intent(in) :: ioInterval + ! Local + real(prec) :: targetTime,tNext + integer :: i,nIO + + this%dt = dt + targetTime = tn + + nIO = int((targetTime-this%t)/ioInterval) + do i = 1,nIO + tNext = this%t+ioInterval + call this%timeIntegrator(tNext) + this%t = tNext + + call this%CalculateEntropy() + call this%ReportEntropy() + call this%ReportMetrics() + + call this%WriteModel() + if(this%tecplot_enabled) then + call this%WriteTecplot() + endif + call this%IncrementIOCounter() + + enddo + + endsubroutine ForwardStep_Model + + subroutine Euler_timeIntegrator(this,tn) + implicit none + class(Model),intent(inout) :: this + real(prec),intent(in) :: tn + ! Local + real(prec) :: tRemain + real(prec) :: dtLim + + dtLim = this%dt ! Get the max time step size from the dt attribute + do while(this%t < tn) + + tRemain = tn-this%t + this%dt = min(dtLim,tRemain) + call this%CalculateTendency() + call this%UpdateSolution() + this%t = this%t+this%dt + + enddo + + this%dt = dtLim + + endsubroutine Euler_timeIntegrator + + subroutine LowStorageRK2_timeIntegrator(this,tn) + implicit none + class(Model),intent(inout) :: this + real(prec),intent(in) :: tn + ! Local + integer :: m + real(prec) :: tRemain + real(prec) :: dtLim + real(prec) :: t0 + + dtLim = this%dt ! Get the max time step size from the dt attribute + do while(this%t < tn) + + t0 = this%t + tRemain = tn-this%t + this%dt = min(dtLim,tRemain) + do m = 1,2 + call this%CalculateTendency() + call this%UpdateGRK2(m) + this%t = t0+rk2_b(m)*this%dt + enddo + + this%t = t0+this%dt + + enddo + + this%dt = dtLim + + endsubroutine LowStorageRK2_timeIntegrator + + subroutine LowStorageRK3_timeIntegrator(this,tn) + implicit none + class(Model),intent(inout) :: this + real(prec),intent(in) :: tn + ! Local + integer :: m + real(prec) :: tRemain + real(prec) :: dtLim + real(prec) :: t0 + + dtLim = this%dt ! Get the max time step size from the dt attribute + do while(this%t < tn) + + t0 = this%t + tRemain = tn-this%t + this%dt = min(dtLim,tRemain) + do m = 1,3 + call this%CalculateTendency() + call this%UpdateGRK3(m) + this%t = t0+rk3_b(m)*this%dt + enddo + + this%t = t0+this%dt - dtLim = this%dt ! Get the max time step size from the dt attribute - do while (this%t < tn) - - tRemain = tn - this%t - this%dt = min(dtLim, tRemain) - call this%CalculateTendency() - call this%UpdateSolution() - this%t = this%t + this%dt - - end do - - this%dt = dtLim - - end subroutine Euler_timeIntegrator - - subroutine LowStorageRK2_timeIntegrator(this, tn) - implicit none - class(Model), intent(inout) :: this - real(prec), intent(in) :: tn - ! Local - integer :: m - real(prec) :: tRemain - real(prec) :: dtLim - real(prec) :: t0 - - dtLim = this%dt ! Get the max time step size from the dt attribute - do while (this%t < tn) - - t0 = this%t - tRemain = tn - this%t - this%dt = min(dtLim, tRemain) - do m = 1, 2 - call this%CalculateTendency() - call this%UpdateGRK2(m) - this%t = t0 + rk2_b(m)*this%dt - end do - - this%t = t0 + this%dt - - end do - - this%dt = dtLim - - end subroutine LowStorageRK2_timeIntegrator - - subroutine LowStorageRK3_timeIntegrator(this, tn) - implicit none - class(Model), intent(inout) :: this - real(prec), intent(in) :: tn - ! Local - integer :: m - real(prec) :: tRemain - real(prec) :: dtLim - real(prec) :: t0 - - dtLim = this%dt ! Get the max time step size from the dt attribute - do while (this%t < tn) - - t0 = this%t - tRemain = tn - this%t - this%dt = min(dtLim, tRemain) - do m = 1, 3 - call this%CalculateTendency() - call this%UpdateGRK3(m) - this%t = t0 + rk3_b(m)*this%dt - end do - - this%t = t0 + this%dt - - end do - - this%dt = dtLim - - end subroutine LowStorageRK3_timeIntegrator - - subroutine LowStorageRK4_timeIntegrator(this, tn) - implicit none - class(Model), intent(inout) :: this - real(prec), intent(in) :: tn - ! Local - integer :: m - real(prec) :: tRemain - real(prec) :: dtLim - real(prec) :: t0 + enddo + + this%dt = dtLim + + endsubroutine LowStorageRK3_timeIntegrator + + subroutine LowStorageRK4_timeIntegrator(this,tn) + implicit none + class(Model),intent(inout) :: this + real(prec),intent(in) :: tn + ! Local + integer :: m + real(prec) :: tRemain + real(prec) :: dtLim + real(prec) :: t0 - dtLim = this%dt ! Get the max time step size from the dt attribute - do while (this%t < tn) + dtLim = this%dt ! Get the max time step size from the dt attribute + do while(this%t < tn) - t0 = this%t - tRemain = tn - this%t - this%dt = min(dtLim, tRemain) - do m = 1, 5 - call this%CalculateTendency() - call this%UpdateGRK4(m) - this%t = t0 + rk4_b(m)*this%dt - end do + t0 = this%t + tRemain = tn-this%t + this%dt = min(dtLim,tRemain) + do m = 1,5 + call this%CalculateTendency() + call this%UpdateGRK4(m) + this%t = t0+rk4_b(m)*this%dt + enddo - this%t = t0 + this%dt + this%t = t0+this%dt - end do + enddo - this%dt = dtLim + this%dt = dtLim - end subroutine LowStorageRK4_timeIntegrator + endsubroutine LowStorageRK4_timeIntegrator -end module SELF_Model +endmodule SELF_Model From 11184a30b7529f4f094d6cecfb7942bbb301917a Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 01:52:51 +0000 Subject: [PATCH 06/11] Add updatedevice after updating sideinfo This resolves #76 --- src/SELF_Mesh_2D_t.f90 | 2 ++ src/SELF_Mesh_3D_t.f90 | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/SELF_Mesh_2D_t.f90 b/src/SELF_Mesh_2D_t.f90 index 31d020d6e..b1ea9c4a2 100644 --- a/src/SELF_Mesh_2D_t.f90 +++ b/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) diff --git a/src/SELF_Mesh_3D_t.f90 b/src/SELF_Mesh_3D_t.f90 index c415c0022..b8d154cd9 100644 --- a/src/SELF_Mesh_3D_t.f90 +++ b/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) From 04d10d5a2eab06c4a70bcd2484c19b62bf084952 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 02:03:29 +0000 Subject: [PATCH 07/11] Use abs(J) in entropy calculation. Resolves #75 . Element-corner node ordering can occur either clockwise or counter-clockwise which can result in a positive or negative jacobian. For the purposes of integration, we want to use a positive values jacobian. --- src/SELF_DGModel2D_t.f90 | 2 +- src/SELF_DGModel3D_t.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SELF_DGModel2D_t.f90 b/src/SELF_DGModel2D_t.f90 index 41a44132b..6045fd216 100644 --- a/src/SELF_DGModel2D_t.f90 +++ b/src/SELF_DGModel2D_t.f90 @@ -344,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 diff --git a/src/SELF_DGModel3D_t.f90 b/src/SELF_DGModel3D_t.f90 index f8c4b13aa..116e34a2c 100644 --- a/src/SELF_DGModel3D_t.f90 +++ b/src/SELF_DGModel3D_t.f90 @@ -342,7 +342,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 From c60141ba5a2bc6dc97876af6aa929aafd92ae86f Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 02:10:40 +0000 Subject: [PATCH 08/11] Fix typo --- src/SELF_DGModel3D_t.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/SELF_DGModel3D_t.f90 b/src/SELF_DGModel3D_t.f90 index 116e34a2c..33e9d77c3 100644 --- a/src/SELF_DGModel3D_t.f90 +++ b/src/SELF_DGModel3D_t.f90 @@ -177,7 +177,8 @@ subroutine ReportMetrics_DGModel3D_t(this) 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 + minv//" , "//maxv + write(output_unit,'(A)',ADVANCE='yes') str enddo call this%ReportUserMetrics() From b248c8fec4496a21a55328175b7e90e3fe23168b Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 03:06:37 +0000 Subject: [PATCH 09/11] Check for finite final entropy --- examples/linear_euler2d_spherical_soundwave_closeddomain.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 b/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 index f2cf58877..9fe9c78e9 100644 --- a/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 +++ b/examples/linear_euler2d_spherical_soundwave_closeddomain.f90 @@ -93,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 From 3d479f632c41456c140f329ab01214fbf660a15c Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 03:08:52 +0000 Subject: [PATCH 10/11] Change "pressure" to "P" --- src/SELF_LinearEuler2D_t.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SELF_LinearEuler2D_t.f90 b/src/SELF_LinearEuler2D_t.f90 index 3d21835d8..755721995 100644 --- a/src/SELF_LinearEuler2D_t.f90 +++ b/src/SELF_LinearEuler2D_t.f90 @@ -101,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 From 3a9a0534955d95537bde664d4adae42895e51dda Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 27 Nov 2024 03:09:14 +0000 Subject: [PATCH 11/11] Use abs(J) in calculateEntropy Resolves #75 on GPUs --- src/gpu/SELF_DGModel2D.f90 | 2 +- src/gpu/SELF_DGModel3D.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/gpu/SELF_DGModel2D.f90 b/src/gpu/SELF_DGModel2D.f90 index 30a651b37..0d5c00042 100644 --- a/src/gpu/SELF_DGModel2D.f90 +++ b/src/gpu/SELF_DGModel2D.f90 @@ -164,7 +164,7 @@ subroutine CalculateEntropy_DGModel2D(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 diff --git a/src/gpu/SELF_DGModel3D.f90 b/src/gpu/SELF_DGModel3D.f90 index 294149cad..fe4e98e07 100644 --- a/src/gpu/SELF_DGModel3D.f90 +++ b/src/gpu/SELF_DGModel3D.f90 @@ -169,7 +169,7 @@ subroutine CalculateEntropy_DGModel3D(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