| | 49 | ==== Here is the bare minimum problem module ==== |
| | 50 | {{{ |
| | 51 | MODULE Problem |
| | 52 | IMPLICIT NONE |
| | 53 | SAVE |
| | 54 | PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, & |
| | 55 | ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep |
| | 56 | PRIVATE |
| | 57 | |
| | 58 | CONTAINS |
| | 59 | |
| | 60 | SUBROUTINE ProblemModuleInit() |
| | 61 | END SUBROUTINE ProblemModuleInit |
| | 62 | |
| | 63 | SUBROUTINE ProblemGridInit(Info) |
| | 64 | TYPE(InfoDef) :: Info |
| | 65 | END SUBROUTINE ProblemGridInit |
| | 66 | |
| | 67 | SUBROUTINE ProblemBeforeStep(Info) |
| | 68 | TYPE(InfoDef) :: Info |
| | 69 | END SUBROUTINE ProblemBeforeStep |
| | 70 | |
| | 71 | SUBROUTINE ProblemAfterStep(Info) |
| | 72 | TYPE(InfoDef) :: Info |
| | 73 | END SUBROUTINE ProblemAfterStep |
| | 74 | |
| | 75 | SUBROUTINE ProblemSetErrFlag(Info) |
| | 76 | TYPE(InfoDef) :: Info |
| | 77 | END SUBROUTINE ProblemSetErrFlag |
| | 78 | |
| | 79 | SUBROUTINE ProblemBeforeGlobalStep(n) |
| | 80 | INTEGER :: n |
| | 81 | END SUBROUTINE ProblemBeforeGlobalStep |
| | 82 | |
| | 83 | END MODULE Problem |
| | 84 | |
| | 85 | |
| | 86 | }}} |
| | 87 | [[CollapsibleEnd]] |
| 201 | | 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 and an overdense spherical clump centered at the origin with radius 1. |
| 202 | | |
| 203 | | Note that during the !ProblemGridInit routine, ghost zones do not need to be initialized (rmbc = 0) - however - during beforestep calculations they should be (if needed). |
| | 240 | 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. |
| 230 | | }}} |
| 231 | | |
| 232 | | Now 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.. |
| | 267 | |
| | 268 | SUBROUTINE ProblemBeforeStep(Info) |
| | 269 | TYPE(InfoDef), POINTER :: Info |
| | 270 | INTEGER :: i,j,k, mbc(3) |
| | 271 | REAL(KIND=xPREC) :: pos(3) |
| | 272 | |
| | 273 | !determine number of ghost zones for each dimension |
| | 274 | mbc=levels(Info%level)%gmbc(levels(Info%level)%step)*merge((/1,1,1/),(/0,0,0/),nDim>=(/1,2,3/)) |
| | 275 | |
| | 276 | ! Initialize wind in leftmost boundary |
| | 277 | DO i=1-mbc(1), Info%mX(1)+mbc(1) |
| | 278 | DO j=1-mbc(2), Info%mX(2)+mbc(2) |
| | 279 | DO k=1-mbc(3), Info%mX(3)+mbc(3) |
| | 280 | pos=CellPos(Info, i, j, k) |
| | 281 | IF (pos(1) < GxBounds(1,1)) THEN |
| | 282 | Info%q(i,j,k,irho) = 1d0 |
| | 283 | Info%q(i,j,k,ivx)=10d0 |
| | 284 | IF (ivy /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivy)=0d0 |
| | 285 | IF (ivz /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivz)=0d0 |
| | 286 | IF (iE /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), iE)=gamma7+50d0 |
| | 287 | END IF |
| | 288 | END DO |
| | 289 | END DO |
| | 290 | END DO |
| | 291 | END SUBROUTINE |
| | 292 | }}} |
| | 293 | |
| | 294 | |
| | 295 | * Note that during the !ProblemGridInit routine, ghost zones do not need to be initialized - but that in !ProblemBeforeStep we need to update any ghost zones. |
| | 296 | * 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. |
| | 297 | |
| | 298 | |
| | 299 | |
| | 300 | 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.. |
| 241 | | PUBLIC ProblemModuleInit, ProblemGridInit, & |
| 242 | | ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep |
| 243 | | PRIVATE REAL(KIND=qPREC) :: rho, radius |
| | 309 | PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, & |
| | 310 | ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep |
| | 311 | PRIVATE |
| | 312 | REAL(KIND=qPREC) :: rho, radius, velocity |
| 278 | | ... |
| 279 | | }}} |
| 280 | | |
| 281 | | 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 astrophysical object and there is a clump object module to assist users in easily creating clumps. For example the above module could be rewritten as: |
| | 347 | |
| | 348 | SUBROUTINE ProblemBeforeStep(Info) |
| | 349 | TYPE(InfoDef), POINTER :: Info |
| | 350 | INTEGER :: i,j,k, mbc(3) |
| | 351 | REAL(KIND=xPREC) :: pos(3) |
| | 352 | |
| | 353 | !determine number of ghost zones for each dimension |
| | 354 | mbc=levels(Info%level)%gmbc(levels(Info%level)%step)*merge((/1,1,1/),(/0,0,0/),nDim>=(/1,2,3/)) |
| | 355 | |
| | 356 | ! Initialize wind in leftmost boundary |
| | 357 | DO i=1-mbc(1), Info%mX(1)+mbc(1) |
| | 358 | DO j=1-mbc(2), Info%mX(2)+mbc(2) |
| | 359 | DO k=1-mbc(3), Info%mX(3)+mbc(3) |
| | 360 | pos=CellPos(Info, i, j, k) |
| | 361 | IF (pos(1) < GxBounds(1,1)) THEN |
| | 362 | Info%q(i,j,k,irho) = 1d0 |
| | 363 | Info%q(i,j,k,ivx)=velocity |
| | 364 | IF (ivy /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivy)=0d0 |
| | 365 | IF (ivz /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), ivz)=0d0 |
| | 366 | IF (iE /= 0) Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3), iE)=gamma7+half*velocity**2 |
| | 367 | END IF |
| | 368 | END DO |
| | 369 | END DO |
| | 370 | END DO |
| | 371 | END SUBROUTINE |
| | 372 | |
| | 373 | SUBROUTINE ProblemAfterStep(Info) |
| | 374 | TYPE(InfoDef) :: Info |
| | 375 | END SUBROUTINE ProblemAfterStep |
| | 376 | |
| | 377 | SUBROUTINE ProblemSetErrFlag(Info) |
| | 378 | TYPE(InfoDef) :: Info |
| | 379 | END SUBROUTINE ProblemSetErrFlag |
| | 380 | |
| | 381 | SUBROUTINE ProblemBeforeGlobalStep(n) |
| | 382 | INTEGER :: n |
| | 383 | END SUBROUTINE ProblemBeforeGlobalStep |
| | 384 | |
| | 385 | END MODULE |
| | 386 | |
| | 387 | }}} |
| | 388 | |
| | 389 | 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 astrophysical object and there is a clump object module to assist users in easily creating clumps. For example the above module could be rewritten using an [AmbientObjects Ambient Object], a [ClumpObjects Clump Object], and a [WindObjects Wind Object]. |
| 291 | | PUBLIC ProblemModuleInit, ProblemGridInit, & |
| 292 | | ProblemBeforeStep, ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep |
| 293 | | PRIVATE REAL(KIND=qPREC) :: rho, radius |
| | 400 | PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, & |
| | 401 | ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep |
| | 402 | PRIVATE |
| | 403 | REAL(KIND=qPREC) :: rho, radius, velocity |
| 315 | | ... |
| 316 | | }}} |
| 317 | | |
| 318 | | So what is going on here? First we've added USE statements for the Clumps and Ambients modules. Then in our ProblemModuleInit() routine we've declared an ambientdef and a clumpdef pointer. Then in CreateAmbient() we create an ambient object with the default properties - density = 1, pressure = 1, velocity = 0, etc... which sets the entire grid to be uniform. Then we create a clump object and modify its properties (density, and radius) before updated the clump object. And we're done. |
| 319 | | |
| 320 | | |
| 321 | | And if we want to get more complicated - we can modify other clump attributes. |
| 322 | | |
| 323 | | |
| 324 | | Each grid (or {{{InfoDef}}} structure) comes with several arrays that describe its geometry: |
| 325 | | |
| 326 | | * '''''mX(3)''''': The grid's overall size. Each element {{{mx(n)}}} represents the number of cells along the ''n''th dimension. |
| 327 | | * '''''xbounds(3,2)''''': The spatial coordinates of the grid's lower and upper bounds in each dimension. These values are ''always'' given in computational units. |
| 328 | | * '''''dX''''': The size of a spatial step. {{{dX}}} can also be thought of as the size of a cell. In AstroBEAR's current form, {{{dX}}} is the same in every dimension and is therefore not actually an array. |
| 329 | | |
| 330 | | To make this example more readable (and therefore easier to debug), we assigned the various array values more Cartesian-sounding names. |
| 331 | | |
| 332 | | Similarly, when calculating the spatial equivalents of an index, we went with {{{x}}}, {{{y}}} and {{{z}}}. The half-step in the position calculations represents the fact that the spatial position is at the ''center'' of the cell, not the left side. |
| 333 | | |
| 334 | | The check on the number of dimensions makes sure that 2D problem does not accidentally get initialized with ghost cells (see below) or a spatial step. This allows us to use the same initialization code for 2D or 3D problems. |
| 335 | | |
| 336 | | The {{{levels()}}} array is a global array structures that contain information specific to each level. You can see in this example that we reference the appropriate level by using the {{{Info%level}}} attribute. |
| | 434 | SUBROUTINE ProblemBeforeStep(Info) |
| | 435 | TYPE(InfoDef) :: Info |
| | 436 | END SUBROUTINE ProblemBeforeStep |
| | 437 | |
| | 438 | SUBROUTINE ProblemAfterStep(Info) |
| | 439 | TYPE(InfoDef) :: Info |
| | 440 | END SUBROUTINE ProblemAfterStep |
| | 441 | |
| | 442 | SUBROUTINE ProblemSetErrFlag(Info) |
| | 443 | TYPE(InfoDef) :: Info |
| | 444 | END SUBROUTINE ProblemSetErrFlag |
| | 445 | |
| | 446 | SUBROUTINE ProblemBeforeGlobalStep(n) |
| | 447 | INTEGER :: n |
| | 448 | END SUBROUTINE ProblemBeforeGlobalStep |
| | 449 | |
| | 450 | END MODULE |
| | 451 | }}} |
| | 452 | |
| | 453 | So what is going on here? First we've added USE statements for the Clumps and Ambients modules. Then in our !ProblemModuleInit() routine we've declared an AmbientDef and a ClumpDef object pointer which we initialize in calls to !CreateAmbient and !CreateClump respectively. We leave the ambient default properties alone (density = 1, pressure = 1, etc...) which sets the entire grid to be at a density of 1 and pressure of 1. Then we create a clump object and modify its properties (density, and radius) before updating the clump object. And we're done. We don't have to worry about dx or mx or cell positions etc... All of that detailed work is done for us. It is important to note that the order that objects are created is the same order they are placed on the grid. So had we created the clump object first - the clump data would have been overwritten by the ambient module and there would be no clump. |
| | 454 | |
| | 455 | And if we want to get more complicated - we can modify other clump attributes. All of the available (and default) options for the clump properties should be documented on the ClumpObjects page. |
| | 456 | |
| | 457 | |
| | 458 | |
