Context Navigation

Changes between Version 51 and Version 52 of ModulesOnAstroBear

Timestamp:: 04/24/17 16:14:12 (8 years ago)
Author:: Jonathan
Comment:: —

Legend:

: Unmodified
: Added
: Removed
: Modified

ModulesOnAstroBear

-              v51
+              v52
 === Using Template Functions ===
 One way to apply our initial/boundary functions to grid patches is to use template subroutines.  This involves writing IC and BC subroutines that adhere to a certain form and then passing them to the ApplyTemplate subroutine.  The ApplyTemplate loops over the grid cells and converts the q-array to primitive form (density, velocity, pressure) before calling the template function (in this case !InitialCondition or !BoundaryCondition.
+The preferred way to apply our initial/boundary functions to grid patches is to use template subroutines.  This involves writing IC and BC subroutines that adhere to a certain form and then passing them to the ApplyTemplate subroutine.  The ApplyTemplate loops over the grid cells and converts the q-array to primitive form (density, velocity, pressure) before calling the template function (in this case !InitialCondition or !BoundaryCondition.
 {{{
 …
 }}}
+A second option is to loop over the cells ourselves
+=== Manually Looping Over Zones ===
 …
 * Also note that all of our checks like (ivy /= 0) make our module "dimension safe" so it can run in 1D, 2D, or 3D and that it could run with an isothermal or ideal EOS.
+Now let's say we want to be able to adjust the density of the clump, the radius of the clump, and the wind velocity at run-time.  To do this we need to declare three variables within our module..
+{{{
+MODULE Problem
+  USE GlobalDeclarations
+  USE PhysicsDeclarations
+  USE DataDeclarations
+  IMPLICIT NONE
+  SAVE
+  PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, &
+         ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
+  PRIVATE
+  REAL(KIND=qPREC) :: rho, radius, velocity
+CONTAINS
+  SUBROUTINE ProblemModuleInit()
+    NAMELIST/ProblemData/ rho, radius, velocity
+    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
+    READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
+    CLOSE(PROBLEM_DATA_HANDLE)
+  END SUBROUTINE
+  SUBROUTINE ProblemGridInit(Info)
+    TYPE(InfoDef) :: Info
+    INTEGER :: i,j,k
+    REAL(KIND=xPREC) :: pos(3)
+    ! Initialize background
+    Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), 1)=1d0
+    Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivx)=0d0
+    IF (ivy /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivy)=0d0
+    IF (ivz /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivz)=0d0
+    IF (iE /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), iE)=gamma7
+    ! Increase density for points inside of clump
+    DO i=1, Info%mX(1)
+      DO j=1, Info%mX(2)
+        DO k=1, Info%mX(3)
+          pos=CellPos(Info, i, j, k)
+          IF (sqrt(sum(pos**2)) < radius) THEN
+            Info%q(i,j,k,irho) = rho
+          END IF
+        END DO
+      END DO
+    END DO
+  END SUBROUTINE
+  SUBROUTINE ProblemBeforeStep(Info)
+    TYPE(InfoDef) :: Info
+    INTEGER :: i,j,k, mbc(3)
+    REAL(KIND=xPREC) :: pos(3)
+    !determine number of ghost zones for each dimension
+    mbc=levels(Info%level)%gmbc(levels(Info%level)%step)*merge((/1,1,1/),(/0,0,0/),nDim>=(/1,2,3/))
+    ! Initialize wind in leftmost boundary
+    DO i=1-mbc(1), Info%mX(1)+mbc(1)
+      DO j=1-mbc(2), Info%mX(2)+mbc(2)
+        DO k=1-mbc(3), Info%mX(3)+mbc(3)
+          pos=CellPos(Info, i, j, k)
+          IF (pos(1) < GxBounds(1,1)) THEN
+            Info%q(i,j,k,irho) = 1d0
+            Info%q(i,j,k,ivx)=velocity
+            IF (iE /= 0) Info%q(i,j,k, iE)=gamma7+half*velocity**2
+            IF (ivy /= 0) Info%q(i,j,k, ivy)=0d0
+            IF (ivz /= 0) Info%q(i,j,k, ivz)=0d0
+          END IF
+        END DO
+      END DO
+    END DO
+  END SUBROUTINE
+  SUBROUTINE ProblemAfterStep(Info)
+    TYPE(InfoDef) :: Info
+  END SUBROUTINE ProblemAfterStep
+  SUBROUTINE ProblemSetErrFlag(Info)
+    TYPE(InfoDef) :: Info
+  END SUBROUTINE ProblemSetErrFlag
+  SUBROUTINE ProblemBeforeGlobalStep(n)
+     INTEGER :: n
+  END SUBROUTINE ProblemBeforeGlobalStep
+END MODULE
+}}}
+=== Using Prebuilt Objects ===
 There are a lot of other ways we could modify this simple example to have the clump be located anywhere, to have a density profile that is smoothed at the edge, to be a different temperature, or move with a particular velocity, etc...  Fortunately, clumps are a commonly used object (as are uniform backgrounds and winds) and there are modules designed to assist users in easily creating clumps, uniform backgrounds, and winds.  For example the above module could be rewritten using an [AmbientObjects Ambient Object], a [ClumpObjects Clump Object], and a [WindObjects Wind Object].
 …
          ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
   PRIVATE
+  REAL(KIND=qPREC) :: rho, radius, velocity
+  REAL(KIND=qPREC) :: rhoAmb, Tamb, rhoClump, RClump, rhoWind, TWind, vWind
+  REAL(KIND=qPREC) :: pWind, pClump, pAmb
 CONTAINS
 …
     CLOSE(PROBLEM_DATA_HANDLE)
+    rhoAmb = rhoAmb / rScale
+    rhoClump = rhoClump / rScale
+    rhoWind = rhoWind / rScale
+    TAmb = TAmb / TempScale
+    TWind = TWind / TempScale
+    RClump = RClump / lScale
+    vWind = vWind / velScale
+    pAmb = rhoAmb * TAmb
+    pWind = rhoWind * TWind
+    pClump = pAmb
     CALL CreateAmbient(Ambient)
+    Ambient%density = rhoAmb
+    Ambient%pressure = pAmb
+    CALL UpdateAmbient(Ambient)
     CALL CreateClump(Clump)
+    Clump%density=rho
+    Clump%radius=radius
+    Clump%density = rhoClump
+    Clump%radius = RClump
+    Clump%pressure = pClump
     CALL UpdateClump(Clump)
     CALL CreateWind(Wind)
+    Wind%velocity=velocity
+    Wind%density = rhoWind
+    Wind%velocity = vWind
+    Wind%temperature = TWind
+    Wind%direction = 1
+    Wind%edge = 1
     CALL UpdateWind(Wind)
 …
 {{{
 &ProblemData
+rho = 10.0
+radius = 1.0
+velocity = 10.0
+rhoAmb = 1.0
+TAmb = 1.0
+rhoClump = 10.0
+rClump = 1.0
+rhoWind = 1.0
+TWind = 1.0
+vWind = 10.0
+/
 }}}