Add gravitational acceleration through environmentals attribute

examples/CompressibleIdealGas2D/HydrostaticAdjustment.f90

@@ -0,0 +1,141 @@
+PROGRAM HydrostaticAdjustment
+
+USE SELF_Constants
+USE SELF_Lagrange
+USE SELF_Mesh
+USE SELF_Geometry
+USE SELF_CLI
+USE SELF_CompressibleIdealGas2D
+
+  IMPLICIT NONE
+
+  INTEGER, PARAMETER :: nvar = 4 ! The number prognostic variables
+
+  REAL(prec) :: dt
+  REAL(prec) :: ioInterval
+  REAL(prec) :: tn
+  INTEGER :: N ! Control Degree
+  INTEGER :: M ! Target degree
+  INTEGER :: quadrature
+  CHARACTER(LEN=self_QuadratureTypeCharLength) :: qChar
+  CHARACTER(LEN=self_QuadratureTypeCharLength) :: integrator
+  LOGICAL :: mpiRequested
+  LOGICAL :: gpuRequested
+
+  REAL(prec) :: referenceEntropy
+  TYPE(Lagrange),TARGET :: interp
+  TYPE(Mesh2D),TARGET :: mesh
+  TYPE(SEMQuad),TARGET :: geometry
+  TYPE(CompressibleIdealGas2D),TARGET :: semModel
+  TYPE(MPILayer),TARGET :: decomp
+  TYPE(CLI) :: args
+  CHARACTER(LEN=255) :: SELF_PREFIX
+  CHARACTER(LEN=500) :: meshfile
+  CHARACTER(LEN=100) :: gravity
+
+
+    CALL get_environment_variable("SELF_PREFIX", SELF_PREFIX)
+    CALL args % Init( TRIM(SELF_PREFIX)//"/etc/cli/default.json")
+    CALL args % LoadFromCLI()
+
+    CALL args % Get_CLI('--output-interval',ioInterval)
+    CALL args % Get_CLI('--end-time',tn)
+    CALL args % Get_CLI('--time-step',dt)
+    CALL args % Get_CLI('--mpi',mpiRequested)
+    CALL args % Get_CLI('--gpu',gpuRequested)
+    CALL args % Get_CLI('--control-degree',N)
+    CALL args % Get_CLI('--control-quadrature',qChar)
+    quadrature = GetIntForChar(qChar)
+    CALL args % Get_CLI('--target-degree',M)
+    CALL args % Get_CLI('--integrator',integrator)
+    CALL args % Get_CLI('--mesh',meshfile)
+
+    IF( TRIM(meshfile) == '')THEN
+      meshfile = TRIM(SELF_PREFIX)//"/etc/mesh/Block2D/Block2D_mesh.h5"
+    ENDIF
+
+    ! Initialize a domain decomposition
+    ! Here MPI is disabled, since scaling is currently
+    ! atrocious with the uniform block mesh
+    CALL decomp % Init(enableMPI=mpiRequested)
+
+    ! Create an interpolant
+    CALL interp % Init(N,quadrature,M,UNIFORM)
+
+    ! Read in mesh file
+    CALL mesh % Read_HOPr(TRIM(meshfile),decomp)
+
+    ! Reset the boundary condition to prescribed
+    CALL mesh % ResetBoundaryConditionType(SELF_BC_NONORMALFLOW)
+
+    ! Generate geometry (metric terms) from the mesh elements
+    CALL geometry % Init(interp,mesh % nElem)
+    CALL geometry % GenerateFromMesh(mesh)
+
+    ! Initialize the semModel
+    CALL semModel % Init(nvar,mesh,geometry,decomp)
+
+    IF( gpuRequested )THEN
+      CALL semModel % EnableGPUAccel()
+      ! Update the device for the whole model
+      ! This ensures that the mesh, geometry, and default state match on the GPU
+      CALL semModel % UpdateDevice()
+    ENDIF
+
+    CALL semModel % SetStaticSTP()
+    CALL semModel % CalculateEntropy()
+    CALL semModel % ReportEntropy()
+    referenceEntropy = semModel % entropy
+
+    ! Uses the initial condition to set the prescribed state
+    ! for model boundary conditions
+    CALL semModel % SetPrescribedSolution()
+
+    CALL semModel % SetGravity( "p = 10.0*y" )
+
+    ! Write the initial condition to file
+    CALL semModel % WriteModel()
+    CALL semModel % WriteTecplot()
+
+    ! Set the time integrator
+    CALL semModel % SetTimeIntegrator(TRIM(integrator))
+
+    ! Set your time step
+    semModel % dt = dt
+
+    !! Forward step the semModel and do the file io
+    CALL semModel % ForwardStep( tn = tn, ioInterval = ioInterval )
+
+    !! Manually write the last semModel state
+    CALL semModel % WriteModel('solution.pickup.h5')
+
+    ! Error checking !
+    IF( semModel % entropy /= semModel % entropy )THEN
+      PRINT*, "Model entropy is not a number"
+      STOP 2
+    ENDIF
+
+    IF( semModel % entropy >= HUGE(1.0_prec) )THEN
+      PRINT*, "Model entropy is infinite."
+      STOP 1
+    ENDIF
+
+    IF( semModel % entropy > referenceEntropy )THEN
+      PRINT*, "Warning : final entropy greater than initial entropy"
+      ! Currently do nothing in this situation, since
+      ! conservative solvers in mapped geometries may
+      ! not be entropy conservative.
+      ! However, throwing this warning will bring some
+      ! visibility
+    ENDIF
+
+    ! Clean up
+    CALL semModel % Free()
+    CALL decomp % Free()
+    CALL geometry % Free()
+    CALL mesh % Free()
+    CALL interp % Free()
+    CALL args % Free()
+    CALL decomp % Finalize()
+
+END PROGRAM HydrostaticAdjustment

src/models/SELF_CompressibleIdealGas2D.f90

@@ -36,6 +36,8 @@ MODULE SELF_CompressibleIdealGas2D
 
     TYPE(MappedScalar2D) :: primitive ! Contains primitive variables 
     TYPE(MappedScalar2D) :: entropyVars ! Contains entropy variables
+    TYPE(MappedScalar2D) :: environmentals ! Functions to describe environmental features (e.g. gravitational potential)
+    TYPE(MappedVector2D) :: environmentalsGradient ! Functions to describe environmental features (e.g. gravitational potential)
 
     TYPE(MappedScalar2D) :: diagnostics
 
@@ -79,6 +81,11 @@ MODULE SELF_CompressibleIdealGas2D
     PROCEDURE,PRIVATE :: SetVelocityFromChar_CompressibleIdealGas2D
     PROCEDURE,PRIVATE :: SetVelocityFromEqn_CompressibleIdealGas2D
 
+    GENERIC :: SetGravity => SetGravityFromChar_CompressibleIdealGas2D,&
+                              SetGravityFromEqn_CompressibleIdealGas2D
+    PROCEDURE,PRIVATE :: SetGravityFromChar_CompressibleIdealGas2D
+    PROCEDURE,PRIVATE :: SetGravityFromEqn_CompressibleIdealGas2D
+
     GENERIC :: SetPrescribedSolution => SetPrescribedSolutionFromChar_CompressibleIdealGas2D,&
                               SetPrescribedSolutionFromEqn_CompressibleIdealGas2D,&
                               SetPrescribedSolutionFromSolution_CompressibleIdealGas2D
@@ -151,12 +158,12 @@ END SUBROUTINE SetBoundaryCondition_CompressibleIdealGas2D_gpu_wrapper
   END INTERFACE
 
   INTERFACE
-    SUBROUTINE Source_CompressibleIdealGas2D_gpu_wrapper(source, solution, N, nVar, nEl) &
+    SUBROUTINE Source_CompressibleIdealGas2D_gpu_wrapper(source, solution, environmentalsGradient, N, nVar, nEl) &
       bind(c,name="Source_CompressibleIdealGas2D_gpu_wrapper")
       USE iso_c_binding
       USE SELF_Constants
       IMPLICIT NONE
-      TYPE(c_ptr) :: source, solution
+      TYPE(c_ptr) :: source, solution, environmentalsGradient
       INTEGER(C_INT),VALUE :: N,nVar,nEl
     END SUBROUTINE Source_CompressibleIdealGas2D_gpu_wrapper
   END INTERFACE
@@ -265,6 +272,8 @@ SUBROUTINE Init_CompressibleIdealGas2D(this,nvar,mesh,geometry,decomp)
     CALL this % prescribedSolution % Init(geometry % x % interp,nvarloc,this % mesh % nElem)
     CALL this % diagnostics % Init(geometry % x % interp,nDiagnostics,this % mesh % nElem)
     CALL this % prescribedDiagnostics % Init(geometry % x % interp,nDiagnostics,this % mesh % nElem)
+    CALL this % environmentals % Init(geometry % x % interp,1,this % mesh % nElem)
+    CALL this % environmentalsGradient % Init(geometry % x % interp,1,this % mesh % nElem)
     CALL this % workSol % Init(geometry % x % interp,nVar,this % mesh % nElem)
     CALL this % prevSol % Init(geometry % x % interp,nVar,this % mesh % nElem)
     CALL this % velocity % Init(geometry % x % interp,1,this % mesh % nElem)
@@ -329,6 +338,10 @@ SUBROUTINE Init_CompressibleIdealGas2D(this,nvar,mesh,geometry,decomp)
     CALL this % primitive % SetUnits(4,"kg/m/s^2")
     CALL this % primitive % SetDescription(4,"Fluid pressure")
 
+    CALL this % environmentals % SetName(1,"gp")
+    CALL this % environmentals % SetUnits(1,"m^2/s^2")
+    CALL this % environmentals % SetDescription(1,"Gravitational Potential")
+
   END SUBROUTINE Init_CompressibleIdealGas2D
 
   SUBROUTINE Free_CompressibleIdealGas2D(this)
@@ -339,6 +352,8 @@ SUBROUTINE Free_CompressibleIdealGas2D(this)
     CALL this % prescribedSolution % Free()
     CALL this % diagnostics % Free()
     CALL this % prescribedDiagnostics % Free()
+    CALL this % environmentals % Free()
+    CALL this % environmentalsGradient % Free()
     CALL this % workSol % Free()
     CALL this % velocity % Free()
     CALL this % compVelocity % Free()
@@ -533,6 +548,69 @@ SUBROUTINE SetVelocityFromChar_CompressibleIdealGas2D(this, eqnChar)
 
   END SUBROUTINE SetVelocityFromChar_CompressibleIdealGas2D
 
+  SUBROUTINE SetGravityFromEqn_CompressibleIdealGas2D(this, eqn)
+  !! Sets the fluid velocity field using the provided equation parser
+  !! The momentum is then updated using the current fluid density field.
+  !! From here, the PreTendency method is called to set other diagnostics
+  !!
+  !! The total energy field is calculated using the internal energy (diagnosed from the
+  !! in-situ temperature) and the new kinetic energy field. 
+    IMPLICIT NONE
+    CLASS(CompressibleIdealGas2D),INTENT(inout) :: this
+    TYPE(EquationParser),INTENT(in) :: eqn
+    ! Local
+    INTEGER :: i,j,iEl,iVar
+
+      CALL this % environmentals % SetEquation(1, eqn % equation)
+
+      CALL this % environmentals % SetInteriorFromEquation( this % geometry, this % t )
+      CALL this % environmentals % BoundaryInterp( gpuAccel = .FALSE. )
+
+      CALL this % environmentals % GradientBR( this % geometry, &
+                                               this % environmentalsGradient, &
+                                               .FALSE. )
+
+      IF( this % gpuAccel )THEN
+        CALL this % environmentals % UpdateDevice()
+        CALL this % environmentalsGradient % UpdateDevice()
+      ENDIF
+
+      CALL this % environmentalsGradient % BoundaryInterp( gpuAccel = this % gpuAccel )
+
+  END SUBROUTINE SetGravityFromEqn_CompressibleIdealGas2D
+
+  SUBROUTINE SetGravityFromChar_CompressibleIdealGas2D(this, eqnChar)
+  !! Sets the fluid velocity field using the provided equation parser
+  !! The momentum is then updated using the current fluid density field.
+  !! From here, the PreTendency method is called to set other diagnostics
+  !!
+  !! The total energy field is calculated using the internal energy (diagnosed from the
+  !! in-situ temperature) and the new kinetic energy field. 
+    IMPLICIT NONE
+    CLASS(CompressibleIdealGas2D),INTENT(inout) :: this
+    CHARACTER(LEN=*),INTENT(in) :: eqnChar
+    ! Local
+    INTEGER :: i,j,iEl,iVar
+    REAL(prec) :: rho, u, v, temperature, internalE, KE
+
+      CALL this % environmentals % SetEquation(1, eqnChar)
+
+      CALL this % environmentals % SetInteriorFromEquation( this % geometry, this % t )
+      CALL this % environmentals % BoundaryInterp( gpuAccel = .FALSE. )
+
+      CALL this % environmentals % GradientBR( this % geometry, &
+                                               this % environmentalsGradient, &
+                                               .FALSE. )
+
+      IF( this % gpuAccel )THEN
+        CALL this % environmentals % UpdateDevice()
+        CALL this % environmentalsGradient % UpdateDevice()
+      ENDIF
+
+      CALL this % environmentalsGradient % BoundaryInterp( gpuAccel = this % gpuAccel )
+
+  END SUBROUTINE SetGravityFromChar_CompressibleIdealGas2D
+
   SUBROUTINE SetPrescribedSolutionFromEqn_CompressibleIdealGas2D(this, eqn) 
     IMPLICIT NONE
     CLASS(CompressibleIdealGas2D),INTENT(inout) :: this
@@ -1080,27 +1158,36 @@ SUBROUTINE Source_CompressibleIdealGas2D(this)
     CLASS(CompressibleIdealGas2D),INTENT(inout) :: this
     ! Local
     INTEGER :: i,j,iEl,iVar
+    REAL(prec) :: rhou, rhov, rho, gx, gy
 
     IF( this % gpuAccel )THEN
 
       CALL Source_CompressibleIdealGas2D_gpu_wrapper( this % source % interior % deviceData, &
               this % solution % interior % deviceData, &
+              this % environmentalsGradient % interior % deviceData, &
               this % source % interp % N, &
               this % source % nVar, &
               this % source % nElem )
 
     ELSE
 
       DO iEl = 1, this % source % nElem
-        DO iVar = 1, this % source % nvar
           DO j = 0, this % source % interp % N
             DO i = 0, this % source % interp % N
 
-              this % source % interior % hostData(i,j,iVar,iEl) = 0.0_prec
+              rhou = this % solution % interior % hostData(i,j,1,iEl)
+              rhov = this % solution % interior % hostData(i,j,2,iEl)
+              rho = this % solution % interior % hostData(i,j,3,iEl)
+              gx = this % environmentalsGradient % interior % hostData(1,i,j,1,iEl)
+              gy = this % environmentalsGradient % interior % hostData(2,i,j,1,iEl) 
+
+              this % source % interior % hostData(i,j,1,iEl) = -rho*gx ! (\rho u)_t = -\rho gx 
+              this % source % interior % hostData(i,j,2,iEl) = -rho*gy ! (\rho v)_t = -\rho gy 
+              this % source % interior % hostData(i,j,3,iEl) = 0.0 ! (\rho )_t = 0 
+              this % source % interior % hostData(i,j,4,iEl) = -rhou*gx-rhou*gy ! (\rho E )_t = -\rho u g_x - \rho u g_
 
             ENDDO
           ENDDO
-        ENDDO
       ENDDO
 
     ENDIF

src/models/hip/SELF_CompressibleIdealGas2D.cpp

@@ -2,23 +2,32 @@
 #include "SELF_HIP_Macros.h"
 #include <cstdio>
 
-__global__ void Source_CompressibleIdealGas2D_gpu(real *source, real *solution, int N, int nVar){
+__global__ void Source_CompressibleIdealGas2D_gpu(real *source, real *solution, real *environmentalsGradient, int N, int nVar){
 
   // Get the array indices from the GPU thread IDs
   size_t iVar = blockIdx.x;
   size_t iEl = blockIdx.y;
   size_t i = threadIdx.x;
   size_t j = threadIdx.y;
 
-    source[SC_2D_INDEX(i,j,iVar,iEl,N,nVar)] = 0.0;
+    real rhou = solution[SC_2D_INDEX(i,j,0,iEl,N,nVar)];
+    real rhov = solution[SC_2D_INDEX(i,j,1,iEl,N,nVar)];
+    real rho = solution[SC_2D_INDEX(i,j,2,iEl,N,nVar)];
+    real gx = environmentalsGradient[VE_2D_INDEX(1,i,j,0,iEl,N,1)];
+    real gy = environmentalsGradient[VE_2D_INDEX(2,i,j,0,iEl,N,1)];
+
+    source[SC_2D_INDEX(i,j,0,iEl,N,nVar)] = -rho*gx;
+    source[SC_2D_INDEX(i,j,1,iEl,N,nVar)] = -rho*gy;
+    source[SC_2D_INDEX(i,j,2,iEl,N,nVar)] = 0.0;
+    source[SC_2D_INDEX(i,j,3,iEl,N,nVar)] = -rhou*gx-rhou*gy;
 
 }
 
 extern "C"
 {
-  void Source_CompressibleIdealGas2D_gpu_wrapper(real **source, real **solution, int N, int nVar, int nEl)
+  void Source_CompressibleIdealGas2D_gpu_wrapper(real **source, real **solution, real **environmentalsGradient, int N, int nVar, int nEl)
   {
-    Source_CompressibleIdealGas2D_gpu<<<dim3(nVar,nEl,1), dim3(N+1,N+1,1), 0, 0>>>(*source, *solution, N, nVar);
+    Source_CompressibleIdealGas2D_gpu<<<dim3(nVar,nEl,1), dim3(N+1,N+1,1), 0, 0>>>(*source, *solution, *environmentalsGradient, N, nVar);
   }
 }