Changes between Version 34 and Version 35 of ModulesOnAstroBear


Ignore:
Timestamp:
01/14/13 21:26:42 (12 years ago)
Author:
Jonathan
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • ModulesOnAstroBear

    v34 v35  
    201201Initializing a grid  involves taking a spatially-constructed problem setup and discretizing it so that it fits nicely in an array.  This process is easiest to explain by dissecting an example, such as the one below, where we are trying to initialize the grid with a uniform background and an overdense spherical clump centered at the origin with radius 1.
    202202
    203 Note that during the !ProblemGridInit routine, ghost zones do not need to be initialized (rmbc = 0) - however - during beforestep calculations they should be.
    204 
    205 {{{
    206 
    207     rmbc=levels(Info%level)%gmbc(levels(Info%level)%step)
    208 
    209     mx=Info%mX(1); dx=levels(Info%level)%dX; xlower=Info%xbounds(1,1)
    210     my=Info%mX(2);                           ylower=Info%xbounds(2,1)
    211     mz=Info%mX(3);                           zlower=Info%xbounds(3,1)
    212 
    213     SELECT CASE(nDim)
    214     CASE(2)
    215        zrmbc=0;mz=1;zl=0;dz=0
    216     CASE(3)
    217        zrmbc=rmbc
    218     END SELECT
    219 
    220     ! Initialize environment
    221     DO k=1-zrmbc, mz+zrmbc
    222         DO j=1-rmbc, my+rmbc
    223             DO i=1-rmbc, mx+rmbc
    224 
    225                 x=(xlower + (REAL(i, xPrec)-half) * dx)
    226                 y=(ylower + (REAL(j, xPrec)-half) * dx)
    227                 z=(zlower + (REAL(k, xPrec)-half) * dx)
    228 
    229                 Info%q(i,j,k,irho) = rho(x,y,z)
    230             END DO
     203Note that during the !ProblemGridInit routine, ghost zones do not need to be initialized (rmbc = 0) - however - during beforestep calculations they should be (if needed).
     204
     205{{{
     206  SUBROUTINE ProblemGridInit(Info)
     207    TYPE(InfoDef), POINTER :: Info
     208    INTEGER :: i,j,k
     209    REAL(KIND=xPREC) :: pos(3)
     210
     211    ! Initialize background
     212    Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), 1)=1d0
     213    Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivx)=0d0
     214    IF (ivy /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivy)=0d0
     215    IF (ivz /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivz)=0d0
     216    IF (iE /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), iE)=gamma7   
     217
     218    ! Increase density for points inside of clump
     219    DO i=1, Info%mX(1)
     220      DO j=1, Info%mX(2)
     221        DO k=1, Info%mX(3)
     222          pos=CellPos(Info, i, j, k)
     223          IF (sqrt(sum(pos**2)) < 1d0) THEN
     224            Info%q(i,j,k,irho) = 10d0
     225          END IF                 
    231226        END DO
     227      END DO
    232228    END DO
    233 }}}
     229  END SUBROUTINE
     230}}}
     231
     232Now let's say we want to be able to adjust the density of the clump and the radius of the clump at run-time.  To do this we need to declare two variables within our module..
     233
     234{{{
     235MODULE Problem
     236  USE GlobalDeclarations
     237  USE PhysicsDeclarations
     238  USE DataDeclarations
     239  IMPLICIT NONE
     240  SAVE
     241  PUBLIC ProblemModuleInit, ProblemGridInit, &
     242       ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
     243  PRIVATE REAL(KIND=qPREC) :: rho, radius
     244
     245CONTAINS
     246
     247  SUBROUTINE ProblemModuleInit()
     248    NAMELIST/ProblemData/ rho, radius
     249    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
     250    READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
     251    CLOSE(PROBLEM_DATA_HANDLE)
     252  END SUBROUTINE
     253 
     254  SUBROUTINE ProblemGridInit(Info)
     255    TYPE(InfoDef), POINTER :: Info
     256    INTEGER :: i,j,k
     257    REAL(KIND=xPREC) :: pos(3)
     258
     259    ! Initialize background
     260    Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), 1)=1d0
     261    Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivx)=0d0
     262    IF (ivy /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivy)=0d0
     263    IF (ivz /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivz)=0d0
     264    IF (iE /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), iE)=gamma7   
     265
     266    ! Increase density for points inside of clump
     267    DO i=1, Info%mX(1)
     268      DO j=1, Info%mX(2)
     269        DO k=1, Info%mX(3)
     270          pos=CellPos(Info, i, j, k)
     271          IF (sqrt(sum(pos**2)) < radius) THEN
     272            Info%q(i,j,k,irho) = rho
     273          END IF                 
     274        END DO
     275      END DO
     276    END DO
     277  END SUBROUTINE
     278}}}
     279
     280
     281
    234282
    235283Each grid (or {{{InfoDef}}} structure) comes with several arrays that describe its geometry: