Blog: BE update: problem.f90

!============================================================================
! This problem places a bonner ebert sphere in the center of the grid. 
!============================================================================


MODULE PROBLEM

  USE DataDeclarations
  USE GlobalDeclarations
  USE PhysicsDeclarations
  USE BE_MODULE
  USE WINDS

  IMPLICIT NONE    ! It's safer to require explicit declarations
  SAVE             ! Save module information
  PRIVATE

  PUBLIC ProblemModuleInit, ProblemGridInit, &
       ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag
  TYPE(pWindDef), DIMENSION(:), ALLOCATABLE :: MyWinds
  INTEGER :: nWinds
  REAL(KIND=qPrec) :: central_rho, xi, clump_rad, rho_out, m_star, at, rho_c, gravo_thermal, time_f, iso_t
  NAMELIST /ProblemData/ central_rho, xi, clump_rad, gravo_thermal

CONTAINS

  SUBROUTINE ProblemModuleInit 
    INTEGER :: i, edge

    ! Read in problem params and vars

    OPEN(UNIT = PROBLEM_DATA_HANDLE, FILE = 'problem.data')
    READ(PROBLEM_DATA_HANDLE, NML=ProblemData)
    CLOSE(PROBLEM_DATA_HANDLE)
    central_rho = central_rho*hmass*Xmu ! now is g/cc
    CALL CalcAmbientParams(central_rho, REAL(xi), clump_rad, M_star, at, rho_out, iso_t)
    gravo_thermal = Sqrt((4*pi*G*central_rho)/(at**2))    ! Gravo_thermal radius -- constant that scales r to xi
    IF (MPI_ID == 0) THEN 

       WRITE (*,*) 'The mass (in solar masses) of your star is:' , M_star/msolar 
       WRITE (*,*) 'The isothermal sound speed (m/s) of your star and surrounding medium is:' , aT
       WRITE (*,*) 'The isothermal temperature (K) of your star and surrounding medium is:', iso_T

       IF ((central_rho/rho_out) < 14.1) THEN  
          WRITE (*,*) 'Your BE sphere IS stable. Rho_c/Rho_o = ', central_rho/rho_out, "(less than 14.1)"
       ELSE 
          WRITE(*,*) '!!Your BE sphere is NOT stable!! Rho_c/Rho_o =', central_rho/rho_out, "(greater than 14.1)"
       END IF
       time_f = clump_rad/at !sound crossing time in cgs
       time_f = time_f/timeScale !sound crossing time in cu
       WRITE(*,*) 'The sound crossing time for your clump in computational units is:', time_f
       WRITE(*,*) 'The central density of your clump (in CGS) is:', central_rho

    END IF

    ALLOCATE(MyWinds(6))
    nWinds = 0
    Do i=1, nDim
       Do edge= 1,2
          IF (Gmthbc(i,edge) ==1) THEN
             nWinds=nWinds +1
             ALLOCATE(MyWinds(nWinds)%p)
             MyWinds(nWinds)%p%dir=i
             MyWinds(nWinds)%p%edge=edge
             CALL InitializeWind(MyWinds(nWinds)%p)
          END IF
       END DO
    END DO

  END SUBROUTINE ProblemModuleInit

  SUBROUTINE ProblemGridInit(Info) !Specific for each info structure
    Type (InfoDef) :: Info
    INTEGER        :: i, j, k, mx, my, mz
    INTEGER        :: mbc, rmbc, zrmbc
    INTEGER        :: nx, ny, nz
    INTEGER        :: ndim
    REAL(KIND=xPrec):: x, y, z, xlower, ylower, zlower,dx, dy, dz, r, gravo_thermal, at


    ! Assigning shorthand notation for variables already defined in Info  
    rmbc=levels(Info%level)%gmbc(levels(Info%level)%step)! Sets up ghost cells and boundaries
    mx=Info%mx(1);  xlower=Info%Xbounds(1,1) 
    my=Info%mx(2);  ylower=Info%Xbounds(2,1)
    mz=Info%mx(3);  zlower=Info%Xbounds(3,1)

    dx=levels(Info%level)%dx

    ! These seem to be "nested" variable declarations a=b=c...?


    SELECT CASE(nDim)
    CASE(2)
       zrmbc=0
       mz=1
       zlower=0
       dz=0  
    CASE(3)
       zrmbc=rmbc
    END SELECT

    ! Initializing the environment
    Info%q(:,:,:,:)=0.d0                                  ! Sanity check -- setting grid to zero
    at = (2*clump_rad*SQRT(PI*G*central_rho))/xi
    gravo_thermal = Sqrt((4*pi*G*central_rho)/(at**2))    ! Gravo_thermal radius -- constant that scales r to xi
    ! Must I allocate an array anywhere for Astrobear 2.0?
    ! i,j,k are indices of cells, and x,y,z are spatial cartesian coords
    Do k=1-zrmbc, mz+zrmbc
       Do j=1-rmbc, my+rmbc
          Do i=1-rmbc, mx+rmbc
             x=(xlower + (REAL(i, xPrec)-half)*dx)
             y=(ylower + (REAL(j, xPrec)-half)*dx)
             z=(zlower + (REAL(k, xPrec)-half)*dx) 
             r = Sqrt(x**2+y**2+z**2) 
             r = (r*lscale)*gravo_thermal                 ! takes r to xi 
             IF (r <= xi) THEN                            ! if r<pc
                Info%q(i,j,k,1) = (central_rho*BE_RHO(REAL(r)))/rscale 
             ELSE 
                Info%q(i,j,k,1) = rho_out/rscale                                              
             END IF
             ! M_star = 8*PI*central_rho*clump_rad**3 *((1 + alpha)/(6 + xi**2) - (6*alpha)/(6 + xi**2)**2 - alpha/(3**(1/alpha)*12*e + xi**2) + D*(2**(-alpha)-1)/ (1+2**(-alpha) * D* xi**2)*(1+D*xi**2))
             Info%q(i,j,k,iE) = (((Info%q(i,j,k,1)*rscale*(aT**2))/gamma)*gamma7)/pScale
          END DO
       END DO
    END DO


    IF ((central_rho/rho_out) < 14.1) THEN 
       WRITE (191,*) '#Your BE sphere IS stable.'
    ELSE 
       WRITE(191,*) '#Your BE sphere is NOT stable!!'
    END IF


  END SUBROUTINE ProblemGridInit

  ! Applies Boundary Conditions
  SUBROUTINE ProblemBeforeStep(INFO)

    TYPE(INFODEF) :: INFO
    INTEGER :: i

  END SUBROUTINE ProblemBeforeStep

  SUBROUTINE ProblemAfterStep(INFO)
    TYPE(INFODEF) :: INFO
    ! No special after step instructions needed; this is a stub.

  END SUBROUTINE ProblemAfterStep

  SUBROUTINE ProblemSetErrFlag(INFO)
    TYPE(INFODEF) :: INFO
    ! No special instructions needed; this is a stub.

  END SUBROUTINE ProblemSetErrFlag


END MODULE PROBLEM