Changes between Version 2 and Version 3 of MultiPhysics

Timestamp:

01/22/13 09:41:20 (12 years ago)

Author:

Jonathan

Comment:

—

MultiPhysics

@@ -19 +19 @@
 ==== MHD ====
 
-Enabling MHD requires two changes to the data files:
- * In {{{global.data}}}, set  the {{{MaintainAuxArrays}}} flag to {{{T}}}.
+To enable MHD:
  * In {{{physics.data}}}, set {{{lMHD}}} to {{{T}}}.
 
-You will also need to initialize the {{{aux}}} arrays in your module file.  In fact, magnetic fields are often initialized by setting the {{{aux}}} values first, and then calculating the cell-centered magnetic fields in {{{q}}} by averaging the {{{aux}}} values on either side of a cell:
+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.  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})]]
+
 
 {{{
-Info%q(i,j,k,iBx) = half * (Info%aux(i,j,k,1) + Info%aux(i+1,j,k,1))
-Info%q(i,j,k,iBy) = half * (Info%aux(i,j,k,1) + Info%aux(i,j+1,k,2))
-Info%q(i,j,k,iBz) = half * (Info%aux(i,j,k,1) + Info%aux(i,j,k+1,3))
+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)
+      x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
+      DO j=1, Info%mX(2)+1
+        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
 }}}