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