| 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 | | |
| 459 | | |
| 460 | | You may have noticed the nested loops don't go from {{{1}}} to {{{mX(n)}}}. This is because the data array size is not the same as the grid size. The size of each grid along dimension ''n'' is {{{mX(n)}}}. The ''data arrays'', however, have a number of "ghost cells" associated with them. The number of ghost cells on the end of a grid is given by the quantity |
| 461 | | |
| 462 | | {{{ |
| 463 | | levels(Info%level)%CoarsenRatio * levels(Info%level)%gmbc |
| 464 | | }}} |
| 465 | | |
| 466 | | These ghost cells only appear along the dimensions of the problem, though, so the data arrays for a 2D problem will not have ghost cells along the third dimension. So the real extents of the data arrays in a 3D problem are: |
| 467 | | |
| 468 | | {{{ |
| 469 | | Info%q(1 - rmbc : mX(1) + rmbc, & |
| 470 | | 1 - rmbc : mX(2) + rmbc, & |
| 471 | | 1 - rmbc : mX(3) + rmbc, & |
| 472 | | NrVars) |
| 473 | | }}} |
| 474 | | |
| 475 | | {{{ |
| 476 | | Info%aux(1 - rmbc : mX(1) + rmbc + 1, & |
| 477 | | 1 - rmbc : mX(2) + rmbc + 1, & |
| 478 | | 1 - rmbc : mX(3) + rmbc + 1, & |
| 479 | | NrVars) |
| 480 | | }}} |
| | 453 | So what is going on here? |
| | 454 | * First we've added USE statements for the Clumps and Ambients modules. |
| | 455 | * Then in our !ProblemModuleInit() routine we've declared an !AmbientDef, a !ClumpDef, and a !WindDef object pointer which we initialize in calls to !CreateAmbient, !CreateClump, and !CreateWind respectively. |
| | 456 | * 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. |
| | 457 | * We modify a few clump properties (density, and radius) before updating the clump object. |
| | 458 | * And finally we set the wind velocity before updating the wind object. |
| | 459 | And we're done. We don't have to worry about dx or mx or cell positions, or the number of ghost zones 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.''' |
| | 460 | |
| | 461 | If we want to get more complicated - we can modify other clump/wind/ambient attributes. All of the available (and default) options for the various objects should be documented on the AstroBearObjects page. |
| | 462 | |
| 486 | | Some modules may need specific regions refined, regardless of whether or not there is any obvious error there. AstroBEAR flags cells for refinement using the array |
| 487 | | |
| 488 | | {{{ |
| 489 | | Info%ErrFlag(1 - rmbc : mX(1) + rmbc, & |
| 490 | | 1 - rmbc : mX(2) + rmbc, & |
| 491 | | 1 - rmbc : mX(2) + rmbc, & |
| 492 | | NrVars) |
| 493 | | }}} |
| 494 | | |
| 495 | | To clear the cell at {{{(i,j,k)}}}, simply set {{{Info%ErrFlag(i,j,k)}}} to 0. An error flag of 0 means that the cell does not ''need'' to be refined. This is in general, not a good idea since a previous routine might have already flagged that cell for refinement with good reason. To mark a cell for refinement, set {{{Info%ErrFlag(i,j,k)}}} to 1. The best place to do this is in the {{{ProblemSetErrFlag()}}} routine; most conventional physical criteria for refinement are already handled by AstroBEAR itself. |
| 496 | | |
| 497 | | See [wiki:ControllingRefinement#ProblemSetErrFlag ProblemSetErrFlag] for more details and an example of how to use this subroutine. |
| 498 | | |
| 499 | | [[BR]] |
| 500 | | |
| 501 | | ==== Sample Module ==== |
| 502 | | |
| 503 | | The best way to understand how a module works is by example. The following example is a module for the [wiki:u/ehansen/RT Rayleigh-Taylor instability]. |
| 504 | | |
| 505 | | [attachment:problem.f90 Sample Problem Module] |
| 506 | | |
| 507 | | [[BR]] |
| | 468 | Some modules may need specific regions refined, regardless of whether or not there is any obvious gradients etc. AstroBEAR flags cells for refinement using the array |
| | 469 | |
| | 470 | {{{ |
| | 471 | Info%ErrFlag(1:mx,1:my,1:mz) |
| | 472 | }}} |
| | 473 | |
| | 474 | To clear the cell at {{{(i,j,k)}}}, simply set {{{Info%ErrFlag(i,j,k)}}} to 0. An error flag of 0 means that the cell does not ''need'' to be refined. This is in general, not a good idea since a previous routine might have already flagged that cell for refinement with good reason. To mark a cell for refinement, set {{{Info%ErrFlag(i,j,k)}}} to 1. The place to do this is in the {{{ProblemSetErrFlag()}}} routine; most conventional physical criteria for refinement are already handled by AstroBEAR itself. For more information see ControllingRefinement for more information as well as [wiki:ControllingRefinement#ProblemSetErrFlag ProblemSetErrFlag] for an example of how to use this subroutine. |
| | 475 | |
| | 476 | [[BR]] |
| | 477 | |
| | 478 | ==== A few notes on magnetic aux fields ==== |
| | 479 | Aux fields are particularly difficult to work with - and if initialized improperly (either form a non-divergenceless physical model or from 2nd order errors due to estimating face averages with face-centered values) will produce a probably small but not insignificant divergence that will stick around for the course of the simulation. The easiest way to avoid divergence in your B-fields is to first calculate the vector potential and then to take the curl discretely. In 2D, this means calculating the value of the vector potential at cell corners (which only has a z-component)- and then differencing them in y to get Bx and differencing them in x to get -By. In 3D, the vector potential should be averaged along each edge of the component parallel to that edge. For example A,,x,, should be averaged along each edge that is parallel to the x-axis. The 2nd order errors due to estimating the average value along an edge by the midpoint will not produce divergence in the resulting B-field. |
| | 480 | |
| | 481 | Let's suppose in we want to initialize the grid with a B-field B,,y,,=sin(x). Well the vector potential would just be A,,z,,=cos(x) and our ProblemGridInit routine would look like: |
| | 482 | {{{ |
| | 483 | SUBROUTINE ProblemGridInit(Info) |
| | 484 | TYPE(InfoDef), POINTER :: Info |
| | 485 | INTEGER :: i,j,k |
| | 486 | REAL(KIND=qPREC) :: pos(3), dx |
| | 487 | dx = levels(Info%level)%dx |
| | 488 | IF (MaintainAuxArrays) THEN |
| | 489 | Info%aux(1:Info%mX(1)+1,1:Info%mX(2), 1:Info%mX(3), 1)=0d0 |
| | 490 | IF (nDim == 3) Info%aux(1:Info%mX(1),1:Info%mX(2), 1:Info%mX(3)+1, 3)=0d0 |
| | 491 | DO i=1, Info%mX(1) |
| | 492 | DO j=1, Info%mX(2)+1 |
| | 493 | DO k=1, Info%mX(3) |
| | 494 | pos=CellPos(Info, i, j, k) |
| | 495 | Info%aux(i,j,k,2)=(cos(pos(1)+half*dx)-cos(pos(1)-half*dx))/dx |
| | 496 | END DO |
| | 497 | END DO |
| | 498 | END DO |
| | 499 | DO i=1, Info%mX(1) |
| | 500 | DO j=1, Info%mX(2) |
| | 501 | DO k=1, Info%mX(3) |
| | 502 | Info%q(i,j,k,iBx)=0d0 |
| | 503 | Info%q(i,j,k,iBy)=half*(Info%aux(i,j,k,2)+Info%aux(i,j+1,k,2)) |
| | 504 | Info%q(i,j,k,iBz)=0d0 |
| | 505 | Info%q(i,j,k,iE)=gamma7+half*Info%q(i,j,k,iBy)**2 |
| | 506 | END DO |
| | 507 | END DO |
| | 508 | END DO |
| | 509 | ELSE |
| | 510 | DO i=1, Info%mX(1) |
| | 511 | DO j=1, Info%mX(2) |
| | 512 | DO k=1, Info%mX(3) |
| | 513 | pos=CellPos(Info, i, j, k) |
| | 514 | Info%q(i,j,k,iBx)=0d0 |
| | 515 | Info%q(i,j,k,iBy)=sin((pos(1)) |
| | 516 | Info%q(i,j,k,iBz)=0d0 |
| | 517 | END DO |
| | 518 | END DO |
| | 519 | END DO |
| | 520 | END IF |
| | 521 | END SUBROUTINE |
| | 522 | }}} |
| | 523 | * Note that if {{{MaintainAuxArrays}}} is false then we don't have to worry about aux fields and we can just initialize the q values with volume averaged (or cell-centered). Currently this only happens in 1D. |
| | 524 | * If {{{MaintainAuxArrays}}} is true then we must initialize all of Info%aux surrounding the '''core''' region. So first we set the Bx and Bz aux fields to zero. |
| | 525 | * Next we iterate over each y-face (notice the +1 in the upper j bound) and take the discrete derivative of the vector potential at the adjacent corners at +- half*dx |
| | 526 | * We then go through and update the cell centered values for the B field using the average of the two adjacent faces - and then we update the total energy to include the magnetic contribution. |
