Changes between Version 4 and Version 5 of MultiPhysics

Timestamp:

01/22/13 16:17:51 (13 years ago)

Author:

Jonathan

Comment:

—

MultiPhysics

@@ -45 +45 @@
     !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)
+    DO i=1, Info%mX(1)+1
       x = Info%xBounds(1,1)+(REAL(i)-1d0)*dx(1)
-      DO j=1, Info%mX(2)+1
+      DO j=1, Info%mX(2)
         y = Info%xBounds(2,1)+(REAL(j)-.5d0)*dx(2)
         DO k=1, Info%mX(3)
@@ -105 +105 @@
 }}}
 
-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 accurately integrated along each edge - but a better approach is to instead initialize the vector potential and then discretely difference it to get the face fields.
+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
+
+}}}
+
+[[Image(Stencil.png, width=400)]]
 
 ==== Cooling ====