| 101 | | ==== Initializing a Grid / Updating boundary conditions ==== |
| 102 | | |
| 103 | | Initializing 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 (density = 1, pressure = 1) and an overdense spherical clump centered at the origin with radius 1. We then want to add a constant wind of density 1, pressure 1, and velocity 10 coming from the left boundary. |
| | 101 | == Reading in Problem parameters and converting to computational units == |
| | 102 | |
| | 103 | Initializing a grid involves setting the fluid values at each cell with spatially averaged values of your initial coniditions. 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 (density = '''rhoAmb''', Temperature = '''Tamb''') and an overdense spherical clump centered at the origin with radius '''R''', density '''rhoClump''' in pressure equilibrium with the ambient. We then want to add a constant wind of density '''rhoWind''', Temperature '''TWind''', and velocity '''Vwind''' coming from the left boundary. |
| | 104 | |
| | 105 | |
| | 106 | First, we need to declare our variables and add them to the !ProblemData Namelist. Then in !ProblemModuleInit, we read them in from the problem.data file and convert all of the variables to computational units. Here I will assume that input variables are in cgs units (with Temperature in Kelvin). We are also going to declare additional derived variables to store the calculated wind pressure, clump pressure, and ambient pressure. |
| | 107 | |
| | 108 | {{{ |
| | 109 | MODULE Problem |
| | 110 | USE GlobalDeclarations |
| | 111 | USE PhysicsDeclarations |
| | 112 | USE DataDeclarations |
| | 113 | IMPLICIT NONE |
| | 114 | SAVE |
| | 115 | PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, & |
| | 116 | ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep |
| | 117 | PRIVATE |
| | 118 | REAL(KIND=qPREC) :: rhoAmb, Tamb, rhoClump, RClump, rhoWind, TWind, vWind |
| | 119 | REAL(KIND=qPREC) :: pWind, pClump, pAmb |
| | 120 | |
| | 121 | CONTAINS |
| | 122 | |
| | 123 | SUBROUTINE ProblemModuleInit() |
| | 124 | NAMELIST/ProblemData/ rhoAmb, Tamb, rhoClump, RClump, rhoWind, TWind, vWind |
| | 125 | OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD") |
| | 126 | READ(PROBLEM_DATA_HANDLE,NML=ProblemData) |
| | 127 | CLOSE(PROBLEM_DATA_HANDLE) |
| | 128 | rhoAmb = rhoAmb / rScale |
| | 129 | rhoClump = rhoClump / rScale |
| | 130 | rhoWind = rhoWind / rScale |
| | 131 | TAmb = TAmb / TempScale |
| | 132 | TWind = TWind / TempScale |
| | 133 | RClump = RClump / lScale |
| | 134 | vWind = vWind / velScale |
| | 135 | |
| | 136 | pAmb = rhoAmb * TAmb |
| | 137 | pWind = rhoWind * TWind |
| | 138 | pClump = pAmb |
| | 139 | |
| | 140 | END SUBROUTINE |
| | 141 | }}} |
| | 142 | |
| | 143 | === Using Template Functions === |
| | 144 | 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. |
| | 145 | |
| | 146 | {{{ |
| | 147 | |
| | 148 | SUBROUTINE ProblemGridInit(Info) |
| | 149 | TYPE(InfoDef) :: Info |
| | 150 | CALL ApplyFunction(InitialCondition, Info) |
| | 151 | END SUBROUTINE ProblemGridInit |
| | 152 | |
| | 153 | SUBROUTINE ProblemBeforeStep(Info) |
| | 154 | TYPE(InfoDef) :: Info |
| | 155 | CALL ApplyFunction(BoundaryCondition, Info) |
| | 156 | END SUBROUTINE ProblemBeforeStep |
| | 157 | |
| | 158 | |
| | 159 | SUBROUTINE InitialCondition(pos, t, dt, q, Ghost) |
| | 160 | REAL(KIND=qPREC), DIMENSION(:) :: pos |
| | 161 | REAL(KIND=qPREC), DIMENSION(:) :: q |
| | 162 | REAL(KIND=qPREC) :: t, dt |
| | 163 | LOGICAL :: Ghost |
| | 164 | IF (sqrt(sum(pos**2)) < RClump) THEN |
| | 165 | q(1) = rhoClump |
| | 166 | q(imom(1:nDim))=0 |
| | 167 | q(iE) = pClump |
| | 168 | ELSE |
| | 169 | q(1) = rhoAmb |
| | 170 | q(imom(1:nDim)) = 0 |
| | 171 | q(iE) = pAmb |
| | 172 | END IF |
| | 173 | END SUBROUTINE InitialCondition |
| | 174 | |
| | 175 | |
| | 176 | SUBROUTINE BoundaryCondition(pos, t, dt, q, Ghost) |
| | 177 | REAL(KIND=qPREC), DIMENSION(:) :: pos |
| | 178 | REAL(KIND=qPREC), DIMENSION(:) :: q |
| | 179 | REAL(KIND=qPREC) :: t, dt |
| | 180 | LOGICAL :: Ghost |
| | 181 | IF (pos(1) < GxBounds(1,1)) THEN |
| | 182 | q(1) = rhoWind |
| | 183 | q(imom(1:nDim))=0 |
| | 184 | q(ivx) = vWind |
| | 185 | q(iE) = pWind |
| | 186 | END IF |
| | 187 | END SUBROUTINE BoundaryCondition |
| | 188 | |
| | 189 | |
| | 190 | }}} |
| | 191 | |
| | 192 | A second option is to loop over the cells ourselves |
| | 193 | |
| 111 | | ! Initialize background |
| 112 | | Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), 1)=1d0 |
| 113 | | Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivx)=0d0 |
| 114 | | IF (ivy /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivy)=0d0 |
| 115 | | IF (ivz /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivz)=0d0 |
| 116 | | IF (iE /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), iE)=gamma7 |
| 117 | | |
| 118 | | ! Increase density for points inside of clump |
| | 201 | ! Loop over cells in patch |
| 123 | | IF (sqrt(sum(pos**2)) < 1d0) THEN |
| 124 | | Info%q(i,j,k,irho) = 10d0 |
| | 206 | IF (sqrt(sum(pos**2)) < RClump) THEN |
| | 207 | Info%q(i,j,k,irho) = rhoClump |
| | 208 | Info%q(i,j,k,imom(1:nDim)) = 0 |
| | 209 | Info%q(i,j,k,iE) = gamma7*pClump |
| | 210 | ELSE |
| | 211 | Info%q(i,j,k,irho) = rhoAmb |
| | 212 | Info%q(i,j,k,imom(1:nDim)) = 0 |
| | 213 | Info%q(i,j,k,iE) = gamma7*pAmb |
