Changes between Version 9 and Version 10 of MultiPhysics


Ignore:
Timestamp:
02/22/17 12:32:19 (8 years ago)
Author:
Jonathan
Comment:

Legend:

Unmodified
Added
Removed
Modified
  • MultiPhysics

    v9 v10  
    1111* [SinkParticles Sink Particles]
    1212
    13 ==== Tracers ====
    14 Using tracers is fairly straightforward in AstroBEAR.  In your !ProblemModuleInit routine you can create additional tracers by calling
    15 {{{
    16  CALL AddTracer(iTracer, 'TracerName')
    17 }}}
    18 where {{{iTracer}}} is an integer that corresponds to the slot in {{{Info%q}}} and !TracerName is an optional string description of the tracer that will show up in visit etc...
    19 ===== Tracers in Objects =====
    20 Most objects such as [ClumpObjects Clumps], [DiskObjects Disks], etc... support tracer fields and will properly initialize the data associated with those.  All you have to do is initialize them as follows:
    21 {{{
    22  CALL CreateClump(Clump)
    23  CALL AddTracer(Clump%iTracer, 'MyClumpTracer')
    24  CALL UpdateClump(Clump)
    25 }}}
    26 
    27 ==== MHD ====
    28 
    29 To enable MHD:
    30  * In {{{physics.data}}}, set {{{lMHD}}} to {{{T}}}.
    31 
    32 If the code is using constrained transport, then the logical flag {{{MaintainAuxArrays}}} will be {{{.True.}}}.  In that case you will need to initialize the face averaged magnetic field values in the {{{aux}}} array in your module file.  The following diagram displays the various aux variables and what positions they correspond to - as well as the positions of the EMF as well as the vector potential.
    33 
    34 [[Image(Stencil.png, width=400)]]
    35 
    36 Here is an example for an initial magnetic field that is [[latex($\mathbf {B}=Bx(x,y,z)\mathbf{e_x} + By(x,y,z) \mathbf{e_y} + Bz(x,y,z) \mathbf{e_z}$)]]
    3713
    3814
    39 {{{
    40 dx = (/(merge(1d0,0d0,nDim >= i),i=1,3)/)*levels(Info%level)%dx
    41 
    42 DO i=1, Info%mX(1)
    43   x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
    44   DO j=1, Info%mX(2)
    45     y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
    46     DO k=1, Info%mX(3)
    47       z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
    48       Info%q(i,j,k,iBx)=Bx(x,y,z)
    49       Info%q(i,j,k,iBy)=By(x,y,z)
    50       Info%q(i,j,k,iBz)=Bz(x,y,z)
    51     END DO
    52   END DO
    53 END DO
    54 
    55 IF (MaintainAuxArrays) THEN
    56   IF (nDim >= 1) THEN
    57     !Note the +1 in the i loop DO statement and the -1d0 in
    58     !the x coordinate calculation instead of -.5d0
    59     DO i=1, Info%mX(1)+1
    60       x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
    61       DO j=1, Info%mX(2)
    62         y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
    63         DO k=1, Info%mX(3)
    64           z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
    65           Info%aux(i,j,k,1)=Bx(x,y,z)
    66         END DO
    67       END DO
    68     END DO
    69 
    70     !Now we can update the cell centered value with the adjacent face average
    71     Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBx)=.5d0 * (&
    72       Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1) + &
    73       Info%aux(2:Info%mX(1)+1,1:Info%mX(2),1:Info%mX(3),1) )
    7415
    7516
    76     IF (nDim >= 2) THEN
    77       !Note the +1 in the j loop DO statement and the -1d0 in
    78       !the y coordinate calculation instead of -.5d0
    79       DO i=1, Info%mX(1)
    80         x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
    81         DO j=1, Info%mX(2)+1
    82           y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
    83           DO k=1, Info%mX(3)
    84             z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
    85             Info%aux(i,j,k,2)=By(x,y,z)
    86           END DO
    87         END DO
    88       END DO
    8917
    90       !Now we can update the cell centered value with the adjacent face average 
    91       Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBy)=.5d0 * (&
    92         Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),2) + &
    93         Info%aux(1:Info%mX(1),2:Info%mX(2)+1,1:Info%mX(3),2) )
    9418
    95       IF (nDim >= 3) THEN
    96        !Note the +1 in the k loop DO statement and the -1d0 in
    97        !the z coordinate calculation instead of -.5d0
    98        DO i=1, Info%mX(1)
    99          x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
    100          DO j=1, Info%mX(2)
    101            y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
    102            DO k=1, Info%mX(3)+1
    103              z = Info%xBounds(3,1)+(REAL(k)-1d0)*dx(3)
    104              Info%aux(i,j,k,3)=Bz(x,y,z)
    105            END DO
    106          END DO
    107        END DO
    10819
    109       !Now we can update the cell centered value with the adjacent face average
    110        Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBz)=.5d0 * (&
    111          Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),3) + &
    112          Info%aux(1:Info%mX(1),1:Info%mX(2),2:Info%mX(3)+1,3) )
    113       END IF
    114     END IF
    115   END IF
    116 END IF
    117 }}}
    118 
    119 Of course the above example works for any field - but does not impose the divergence free constraint. Presumably the functions {{{Bx}}}, {{{By}}}, and {{{Bz}}} would take care of that to the degree they are accurate face averages - but a better approach is to instead approximate  the vector potential at cell edges and then discretely difference it to get the face average fields.
    120 
    121 {{{
    122 IF (nDim == 2) THEN
    123   !In 2D the only component of the vector potential that matters is Az
    124   ALLOCATE(A(1:Info%mX(1)+1,1:Info%mX(2)+1,1,1))
    125   DO i=1, Info%mX(1)+1
    126     x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
    127     DO j=1, Info%mX(2)+1
    128       y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
    129       A(i,j,1,1)=Az(x,y,z)
    130     END DO
    131   END DO
    132 
    133   !Then take the curl of A to get B
    134   FORALL(i=1:Info%mX(1)+1,j=1:Info%mX(2))
    135     Info%aux(i,j,1,1)=+(A(i,j+1,1,1)-A(i,j,1,1))/dx(2)
    136   END FORALL
    137   FORALL(i=1:Info%mX(1),j=1:Info%mX(2)+1)
    138     Info%aux(i,j,1,2)=-(A(i+1,j,1,1)-A(i,j,1,1))/dx(1)
    139   END FORALL
    140 
    141 ELSEIF (nDim == 3) THEN
    142   !Now the vector potential is a 3 component vector
    143   ALLOCATE(A(1:Info%mX(1)+1,1:Info%mX(2)+1,1:Info%mX(3)+1,3))
    144 
    145   !We could do this with three separate loops to avoid calculating
    146   !A(:,:,Info%mX(3)+1,3),A(:,Info%mX(2)+1,:,2), and A(Info%mX(1)+1,:,:,1)
    147   !Since these do not get used in taking the curl...
    148 
    149   DO i=1, Info%mX(1)+1
    150     x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
    151     DO j=1, Info%mX(2)+1
    152       y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
    153       DO k=1, Info%mX(3)+1
    154         z = Info%xBounds(3,1)+(REAL(k)-1d0)*dx(3)
    155         A(i,j,k,1)=Ax(x+half*dx(1),y,z)
    156         A(i,j,k,2)=Ay(x,y+half*dx(2),z)
    157         A(i,j,k,3)=Az(x,y,z+half*dx(3))
    158       END DO
    159     END DO
    160   END DO
    161 
    162   !Then take the curl of A to get B
    163   FORALL(i=1:Info%mX(1)+1,j=1:Info%mX(2),k=1:Info%mX(3))
    164     Info%aux(i,j,k,1)=+(A(i,j+1,k,3)-A(i,j,k,3))/dx(2)-(A(i,j,k+1,2)-A(i,j,k,2))/dx(3)
    165   END FORALL
    166   FORALL(i=1:Info%mX(1),j=1:Info%mX(2)+1,k=1:Info%mX(3))
    167     Info%aux(i,j,k,2)=+(A(i,j,k+1,1)-A(i,j,k,1))/dx(3)-(A(i+1,j,k,3)-A(i,j,k,3))/dx(1)
    168   END FORALL
    169   FORALL(i=1:Info%mX(1),j=1:Info%mX(2),k=1:Info%mX(3)+1)
    170     Info%aux(i,j,k,3)=+(A(i+1,j,k,2)-A(i,j,k,2))/dx(1)-(A(i,j+1,k,1)-A(i,j,k,1))/dx(2)
    171   END FORALL
    172 END IF
    173 
    174 }}}
    175 
    176 ==== Cooling ====
    177 
    178 Two things are required to turn on cooling:  the {{{lCooling}}} flag to indicate that cooling is active in this simulation, and {{{iCooling}}} to specify the type of cooling to use.  These values are usually included in the {{{problem.data}}} file.  The user must also create a cooling object in {{{ProblemModuleInit()}}} to manage the cooling settings.  An example of cooling  object creation can be seen below:
    179 
    180 {{{
    181     IF(iCooling>0) THEN
    182        IF (.NOT. lRestart) THEN
    183            ! see sources/cooling.f90::CreateCoolingObject for
    184            ! default values of a cooling source term
    185            CALL CreateCoolingObject(coolingobj)
    186        ELSE
    187            coolingobj => firstcoolingobj
    188        END IF
    189     END IF
    190 
    191     coolingobj%iCooling=iCooling
    192     SELECT CASE(iCooling) ! cases defined in sources/cooling.f90
    193     CASE(NoCool)
    194     CASE(AnalyticCool)
    195        coolingobj%alpha=alpha
    196        coolingobj%beta=beta
    197     CASE(DMCool)
    198     CASE(IICool)
    199     CASE DEFAULT
    200     END SELECT
    201 
    202     coolingobj%floortemp=1d0
    203     coolingobj%mintemp=0.001
    204 }}}
    205 
    206 The {{{.NOT. lRestart}}} conditional prevents AstroBEAR from creating a new cooling object on restarts; this is because the cooling objects will be read in from the restart files.
    207 
    208 ==== Self-Gravity ====
    209 
    210 AstroBEAR uses the [https://computation.llnl.gov/casc/hypre/software.html hypre] library to solve the self-gravity equations.  To use self-gravity:
    211 
    212 1.  Look for the {{{HYPREFLAG}}} variable in {{{Makefile.inc}}} and make sure that it is set to {{{1}}}.
    213 2.  Set the {{{lSelfGravity}}} flag in your {{{physics.data}}} file and set it to {{{T}}}.
    214 
    215 Hypre will automatically initialize the potential field using the density.  The only caveat is that the initial density cannot be uniform.  When the density is uniform, hypre produces a [http://mathworld.wolfram.com/SingularMatrix.html singular matrix] that it can't solve.  Fortunately, a small density perturbation takes care of this problem without substantially affecting the dynamics of the domain.  AstroBEAR comes with a Perturbation object type that can be used for this.
    216 
    217 ==== Sink Particles ====
    218 
    219 The ability to form [SinkParticles sink particles] in AstroBEAR is tied to self-gravity.  To simply enable sink particles:
    220 
    221 1.  Look for the {{{HYPREFLAG}}} variable in {{{Makefile.inc}}} and make sure that it is set to {{{1}}}.
    222 2.  Set the {{{lSelfGravity}}} flag in your {{{physics.data}}} file and set it to {{{T}}}.
    223 
    224 If you just want your simulation to have the option of forming sink particles, no further action is required.  If you want to start your simulation off with sink particles, then you will have to create one in {{{problem.f90::ProblemModuleInit()}}}:
    225 
    226 {{{
    227       NAMELIST /ProblemData/ nParticles
    228       NAMELIST /ParticleData/ mass,xloc,vel
    229       OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
    230       READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
    231 
    232       IF (.NOT. lRestart) THEN
    233 
    234          DO i=1,nParticles
    235 
    236             READ(PROBLEM_DATA_HANDLE,NML=ParticleData)
    237             NULLIFY(Particle)
    238             CALL CreateParticle(Particle)
    239             Particle%mass=mass
    240             Particle%xloc=xloc
    241             Particle%vel=vel         
    242             CALL AddSinkParticle(Particle)
    243  
    244          END DO
    245 
    246          CLOSE(PROBLEM_DATA_HANDLE)     
    247          OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='restart.data', STATUS="UNKNOWN")
    248          WRITE(PROBLEM_DATA_HANDLE,NML=RestartData)
    249          CLOSE(PROBLEM_DATA_HANDLE)
    250 
    251       END IF
    252 }}}
    253 
    254 Depending on the features of your simulation, more [AstroBearObjects objects] might have to be declared in conjunction with the sink particle.  The {{{.NOT. lRestart}}} conditional is important, as it prevents AstroBEAR from adding the same particle again on a restart.