Additional Physics

AstroBEAR supports hydrodynamics and AMR by default. Other physical processes such as magnetic fields and source terms require extra overhead, so they must be enabled by the user.

• Tracers
• MHD
• Cooling
• Radiation driving
• Point Gravity sources
• Self gravity
• Sink Particles

Tracers

Using tracers is fairly straightforward in AstroBEAR. In your ProblemModuleInit routine you can create additional tracers by calling

 CALL AddTracer(iTracer, 'TracerName')

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…

Tracers in Objects

Most objects such as Clumps, Disks, etc… support tracer fields and will properly initialize the data associated with those. All you have to do is initialize them as follows:

 CALL CreateClump(Clump)
 CALL AddTracer(Clump%iTracer, 'MyClumpTracer')
 CALL UpdateClump(Clump)

MHD

To enable MHD:

• In physics.data, set lMHD to T.

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.

Here is an example for an initial magnetic field that is \mathbf {B}=Bx(x,y,z)\mathbf{e_x} + By(x,y,z) \mathbf{e_y} + Bz(x,y,z) \mathbf{e_z}

dx = (/(merge(1d0,0d0,nDim >= i),i=1,3)/)*levels(Info%level)%dx

DO i=1, Info%mX(1)
  x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
  DO j=1, Info%mX(2)
    y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
    DO k=1, Info%mX(3)
      z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
      Info%q(i,j,k,iBx)=Bx(x,y,z)
      Info%q(i,j,k,iBy)=By(x,y,z)
      Info%q(i,j,k,iBz)=Bz(x,y,z)
    END DO
  END DO
END DO

IF (MaintainAuxArrays) THEN
  IF (nDim >= 1) THEN
    !Note the +1 in the i loop DO statement and the -1d0 in 
    !the x coordinate calculation instead of -.5d0
    DO i=1, Info%mX(1)+1
      x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
      DO j=1, Info%mX(2)
        y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
        DO k=1, Info%mX(3)
          z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
          Info%aux(i,j,k,1)=Bx(x,y,z)
        END DO
      END DO
    END DO

    !Now we can update the cell centered value with the adjacent face average
    Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBx)=.5d0 * (&
      Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),1) + &
      Info%aux(2:Info%mX(1)+1,1:Info%mX(2),1:Info%mX(3),1) )


    IF (nDim >= 2) THEN
      !Note the +1 in the j loop DO statement and the -1d0 in 
      !the y coordinate calculation instead of -.5d0
      DO i=1, Info%mX(1)
        x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
        DO j=1, Info%mX(2)+1
          y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
          DO k=1, Info%mX(3)
            z = Info%xBounds(3,1)+(REAL(k)-.5d0)*dx(3)
            Info%aux(i,j,k,2)=By(x,y,z)
          END DO
        END DO
      END DO

      !Now we can update the cell centered value with the adjacent face average  
      Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBy)=.5d0 * (&
        Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),2) + &
        Info%aux(1:Info%mX(1),2:Info%mX(2)+1,1:Info%mX(3),2) )

      IF (nDim >= 3) THEN
       !Note the +1 in the k loop DO statement and the -1d0 in 
       !the z coordinate calculation instead of -.5d0
       DO i=1, Info%mX(1)
         x = Info%xBounds(1,1)+(REAL(i)-.5d0)*dx(1)
         DO j=1, Info%mX(2)
           y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
           DO k=1, Info%mX(3)+1
             z = Info%xBounds(3,1)+(REAL(k)-1d0)*dx(3)
             Info%aux(i,j,k,3)=Bz(x,y,z)
           END DO
         END DO
       END DO

      !Now we can update the cell centered value with the adjacent face average
       Info%q(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),iBz)=.5d0 * (&
         Info%aux(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3),3) + &
         Info%aux(1:Info%mX(1),1:Info%mX(2),2:Info%mX(3)+1,3) )
      END IF
    END IF
  END IF
END IF

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.

IF (nDim == 2) THEN
  !In 2D the only component of the vector potential that matters is Az
  ALLOCATE(A(1:Info%mX(1)+1,1:Info%mX(2)+1,1,1))
  DO i=1, Info%mX(1)+1
    x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
    DO j=1, Info%mX(2)+1
      y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
      A(i,j,1,1)=Az(x,y,z)
    END DO
  END DO

  !Then take the curl of A to get B
  FORALL(i=1:Info%mX(1)+1,j=1:Info%mX(2))
    Info%aux(i,j,1,1)=+(A(i,j+1,1,1)-A(i,j,1,1))/dx(2)
  END FORALL
  FORALL(i=1:Info%mX(1),j=1:Info%mX(2)+1)
    Info%aux(i,j,1,2)=-(A(i+1,j,1,1)-A(i,j,1,1))/dx(1)
  END FORALL

ELSEIF (nDim == 3) THEN
  !Now the vector potential is a 3 component vector
  ALLOCATE(A(1:Info%mX(1)+1,1:Info%mX(2)+1,1:Info%mX(3)+1,3))

  !We could do this with three separate loops to avoid calculating 
  !A(:,:,Info%mX(3)+1,3),A(:,Info%mX(2)+1,:,2), and A(Info%mX(1)+1,:,:,1)
  !Since these do not get used in taking the curl...

  DO i=1, Info%mX(1)+1
    x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
    DO j=1, Info%mX(2)+1
      y = Info%xBounds(2,1)+(REAL(j)-1d0)*dx(2)
      DO k=1, Info%mX(3)+1
        z = Info%xBounds(3,1)+(REAL(k)-1d0)*dx(3)
        A(i,j,k,1)=Ax(x+half*dx(1),y,z)
        A(i,j,k,2)=Ay(x,y+half*dx(2),z)
        A(i,j,k,3)=Az(x,y,z+half*dx(3))
      END DO
    END DO
  END DO

  !Then take the curl of A to get B
  FORALL(i=1:Info%mX(1)+1,j=1:Info%mX(2),k=1:Info%mX(3))
    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)
  END FORALL
  FORALL(i=1:Info%mX(1),j=1:Info%mX(2)+1,k=1:Info%mX(3))
    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)
  END FORALL
  FORALL(i=1:Info%mX(1),j=1:Info%mX(2),k=1:Info%mX(3)+1)
    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)
  END FORALL
END IF

Cooling

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:

    IF(iCooling>0) THEN
       IF (.NOT. lRestart) THEN
           ! see sources/cooling.f90::CreateCoolingObject for
           ! default values of a cooling source term
           CALL CreateCoolingObject(coolingobj)
       ELSE
           coolingobj => firstcoolingobj
       END IF
    END IF

    coolingobj%iCooling=iCooling
    SELECT CASE(iCooling) ! cases defined in sources/cooling.f90
    CASE(NoCool)
    CASE(AnalyticCool)
       coolingobj%alpha=alpha
       coolingobj%beta=beta
    CASE(DMCool)
    CASE(IICool)
    CASE DEFAULT
    END SELECT

    coolingobj%floortemp=1d0
    coolingobj%mintemp=0.001

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.

Self-Gravity

AstroBEAR uses the ​hypre library to solve the self-gravity equations. To use self-gravity:

1. Look for the HYPREFLAG variable in Makefile.inc and make sure that it is set to 1.
2. Set the lSelfGravity flag in your physics.data file and set it to T.

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 ​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.

Sink Particles

The ability to form sink particles in AstroBEAR is tied to self-gravity. To simply enable sink particles:

1. Look for the HYPREFLAG variable in Makefile.inc and make sure that it is set to 1.
2. Set the lSelfGravity flag in your physics.data file and set it to T.

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():

      NAMELIST /ProblemData/ nParticles
      NAMELIST /ParticleData/ mass,xloc,vel
      OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
      READ(PROBLEM_DATA_HANDLE,NML=ProblemData)

      IF (.NOT. lRestart) THEN

         DO i=1,nParticles

            READ(PROBLEM_DATA_HANDLE,NML=ParticleData)
            NULLIFY(Particle)
            CALL CreateParticle(Particle)
            Particle%mass=mass
            Particle%xloc=xloc
            Particle%vel=vel         
            CALL AddSinkParticle(Particle)
 
         END DO

         CLOSE(PROBLEM_DATA_HANDLE)      
         OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='restart.data', STATUS="UNKNOWN")
         WRITE(PROBLEM_DATA_HANDLE,NML=RestartData)
         CLOSE(PROBLEM_DATA_HANDLE)

      END IF

Depending on the features of your simulation, more 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.