ControllingRefinement: problem.f90

!---------------------------------------------------------------
! This is the problem module for the Rayleigh-Taylor instability.
!---------------------------------------------------------------
MODULE Problem

  ! In general, DataDeclarations is always used since it defines the InfoDef type
  ! PhysicsDeclarations and GlobalDeclarations are almost always used
  ! To see exactly what these "sub-modules" define, look at their .f90 files in the astrobear directory
  USE DataDeclarations
  USE PhysicsDeclarations
  USE GlobalDeclarations
  USE SourceDeclarations
  IMPLICIT NONE
  SAVE
  
  ! This PUBLIC statement must be present
  PUBLIC ProblemModuleInit, ProblemGridInit, &
     ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag
  
  ! Some variables are used in more than one SUBROUTINE so they are defined here
  ! qPREC is a type of precision defined in GlobalDeclarations
  REAL(KIND=qPREC), PUBLIC :: DensityAbove, DensityBelow, Amplitude, atwood_number, wave_number, lambda_analytic

CONTAINS

  ! Initialize module variables
  SUBROUTINE ProblemModuleInit()

    ! Define variables and namelists
    INTEGER :: OpenStatus
    NAMELIST /ProblemData/ DensityAbove, DensityBelow, Amplitude
    
    ! Open problem.data and read in data
    ! Problem_Data_Handle is just an integer defined in GlobalDeclarations
    ! OpenStatus variable could be used to check if file is opening correctly
    OPEN(UNIT=Problem_Data_Handle, FILE='problem.data', STATUS="OLD",IOSTAT=OpenStatus)
    READ(Problem_Data_Handle, NML=ProblemData)
    CLOSE(Problem_Data_Handle)
  END SUBROUTINE ProblemModuleInit

  ! Initialize data arrays
  SUBROUTINE ProblemGridInit(Info)

    ! User defined variables, note which ones must be REAL and which are INTEGER
    TYPE(InfoDef) :: Info
    INTEGER :: rmbc, zrmbc, mx, my, mz, i, j, k
    REAL(KIND=qPREC) :: x, y, z, rho, P, Lx, Ly, Lz, vy
    REAL(KIND=qPREC) :: dx, dz, xlower, ylower, zlower

    ! Initializes the q-array
    Info%q=0
   
    ! levels is defined in GlobalDeclarations
    rmbc=levels(Info%level)%gmbc(levels(Info%level)%step)

    ! Assign some useful values with convenient names for ease of use and readability
    ! Gxbounds is defined in GlobalDeclarations
    mx=Info%mX(1); dx=levels(Info%level)%dX; xlower=Info%xbounds(1,1); Lx=Gxbounds(1,2)-Gxbounds(1,1)
    my=Info%mX(2); dz=levels(Info%level)%dX; ylower=Info%xbounds(2,1); Ly=Gxbounds(2,2)-Gxbounds(2,1)
    mz=Info%mX(3);                           zlower=Info%xbounds(3,1); Lz=Gxbounds(3,2)-Gxbounds(3,1)

    ! Allows for 2D or 3D simulations
    ! Note that in 2D, Lz is explicitly set to 1 so there is no divide by zero in the perturbation
    SELECT CASE(nDim)
    CASE(2)
       zrmbc=0;mz=1;zlower=0;dz=0;Lz=1
    CASE(3)
       zrmbc=rmbc
    END SELECT



    ! Initialize the grid
    DO k=1-zrmbc, mz+zrmbc
        DO j=1-rmbc, my+rmbc
            DO i=1-rmbc, mx+rmbc

                ! Cell-to-space conversion (half is defined in GlobalDeclarations)
                x=(xlower + (REAL(i) - half) * dx)
                y=(ylower + (REAL(j) - half) * dx)
                z=(zlower + (REAL(k) - half) * dz)
                
                ! Define density profile
                IF (y < 0) THEN
                   rho = DensityBelow
                ELSE
                   rho = DensityAbove
                END IF

                ! Define pressure gradient
                P = 2.5 - rho * mCentral * y

                ! Define perturbation (Pi is defined in PhysicsDeclarations)
                vy = Amplitude * (1+COS(2*Pi*x/Lx)) * (1+COS(2*Pi*y/Ly)) * (1+COS(2*Pi*z/Lz)) / 8

                ! Put information into q-array (gamma, ivy, iE are defined in PhysicsDeclarations)
                Info%q(i, j, k, 1) = rho
                Info%q(i, j, k, ivy) = rho * vy
                Info%q(i, j, k, iE) = P / (gamma - 1.0) + half * rho * vy**2
            END DO
        END DO
    END DO

    ! For analysis...calculates the analytic growth rate and saves it in a file
    atwood_number = (DensityAbove - DensityBelow) / (DensityAbove + DensityBelow)
    wave_number = 2*Pi/Lx
    lambda_analytic = SQRT(atwood_number * wave_number * mCentral)
    OPEN(UNIT=82,FILE='GrowthRate.data',STATUS='REPLACE')
    WRITE(82,*) lambda_analytic
    CLOSE(82)

  END SUBROUTINE ProblemGridInit

  ! Place any pre-processing operations here (This is what is meant by leaving a SUBROUTINE as a stub)
  SUBROUTINE ProblemBeforeStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemBeforeStep

  ! Place any post-processing operations here
  SUBROUTINE ProblemAfterStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemAfterStep

  ! Can be used to set additional refinement
  SUBROUTINE ProblemSetErrFlag(Info)
    TYPE(InfoDef) :: Info
    REAL(KIND=qPREC) :: pos
    REAL(KIND=qPREC), PARAMETER :: lambda = 0.6472
    REAL(KIND=qPREC), DIMENSION(0:MaxLevel) :: mybuffer
    REAL(KIND=qPREC), DIMENSION(:,:), POINTER :: offsets
    REAL(KIND=qPREC), DIMENSION(3,2) :: xBounds
    INTEGER, DIMENSION(3,2) :: mS
    INTEGER, DIMENSION(:,:,:), POINTER :: mSs
    INTEGER :: i, nOverlaps
    mybuffer(0:MaxLevel)=(/5d-2,25d-3,125d-4/)
    xBounds=GxBounds
    pos=0.003*EXP(lambda*levels(Info%level)%tnow) + mybuffer(Info%level)
    xBounds(2,:)=(/-pos,+pos/)
    CALL CalcPhysicalOverlaps(Info, xBounds, mSs, nOverlaps, offsets, IEVERYWHERE, lHydroPeriodic,0)
    IF (nOverlaps > 0) THEN
       DO i=1,nOverlaps
          mS=mSs(i,:,:)
          Info%ErrFlag(mS(1,1):mS(1,2), mS(2,1):mS(2,2), mS(3,1):mS(3,2)) = 1
       END DO
       DEALLOCATE(mSs, offsets)
       NULLIFY(mSs, offsets)
    END IF
  END SUBROUTINE ProblemSetErrFlag

END MODULE Problem