tracking down error with potential source calculation...

fixed error in grid size determination in calculus module

more updates to deal with derivatives

src/numerical/calculus/div.f90

@@ -186,7 +186,7 @@
   end do
 end do
 
-if (lx2>1) then   !only if the x2-direction is not null
+if (x%lx2>1) then   !only if the x2-direction is not null, this should be based on the grid data and not the data passed into this function
   do ix3=1,lx3
     do ix1=1,lx1
       div3D_curv_23(ix1,1,ix3)=div3D_curv_23(ix1,1,ix3)+ &
@@ -205,7 +205,7 @@
 end if
 
 
-if (lx3>1) then    !only if non-singleton 3rd dimension
+if (x%lx3>1) then    !only if non-singleton 3rd dimension (I'm not sure this will ever be the case due to dimension swapping)
   do ix2=1,lx2
     do ix1=1,lx1
       div3D_curv_23(ix1,ix2,1)=div3D_curv_23(ix1,ix2,1)+ &
@@ -228,4 +228,4 @@
 
 end procedure div3D_curv_23
 
-end submodule div
+end submodule div

src/numerical/calculus/gradient.f90

@@ -226,7 +226,7 @@
 lx2=size(f,2)
 lx3=size(f,3)
 
-if (lx2>1) then    !if we have a singleton dimension then we are doing a 2D run and the derivatives in this direction are zero
+if (x%lx2>1) then    !if we have a singleton dimension then we are doing a 2D run and the derivatives in this direction are zero
 
   !ERROR CHECKING TO MAKE SURE DIFFRENCING IS DONE OVER A CONSISTENTLY-SIZED GRID
   if (lx1 /= ubnd1-lbnd1+1 .or. lx2 /= ubnd2-lbnd2+1 .or. lx3 /= ubnd3-lbnd3+1) then
@@ -344,7 +344,7 @@
 lx3=size(f,3)
 
 
-if (lx3>1) then    !only differentiate if we have non-singleton dimension, otherwise set to zero
+if (x%lx3>1) then    !only differentiate if we have non-singleton dimension, otherwise set to zero
   !ERROR CHECKING TO MAKE SURE DIFFRENCING IS DONE OVER A CONSISTENTLY-SIZED GRID
   if (lx1 /= ubnd1-lbnd1+1 .or. lx2 /= ubnd2-lbnd2+1 .or. lx3 /= ubnd3-lbnd3+1) then
     error stop '!!!  Inconsistent array and mesh sizes in gradient function.'   !just bail on it and let the user figure it out

src/numerical/potential/potential_comm_mumps.f90

@@ -6,6 +6,8 @@ module potential_comm
 !!
 !! NOTE THAT ONLY THE CURVILINEAR FUNCTION ARE UP TO DATE.
 
+use, intrinsic :: ieee_arithmetic
+
 use phys_consts, only: wp, pi, lsp, debug
 use grid, only: flagswap, gridflag
 use mesh, only: curvmesh
@@ -156,7 +158,6 @@ subroutine electrodynamics_curv(it,t,dt,nn,vn2,vn3,Tn,sourcemlat,ns,Ts,vs1,B1,vs
 call cpu_time(tstart)
 call conductivities(nn,Tn,ns,Ts,vs1,B1,sig0,sigP,sigH,muP,muH,muPvn,muHvn)
 
-
 if (flagcap/=0) then
   call capacitance(ns,B1,flagcap,incap)
   if (it==1) then     !check that we don't have an unsupported grid type for doing capacitance

src/numerical/potential/potential_root.f90

@@ -140,6 +140,11 @@
              J3halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
 srcterm=divtmp(1:lx1,1:lx2,1:lx3)
 if (debug) print *, 'Root has computed background field source terms...',minval(srcterm), maxval(srcterm)
+
+!print*, myid, any(ieee_is_nan(J1halo(0:lx1+1,1:lx2,1:lx3))), &
+!              any(ieee_is_nan(J2halo(1:lx1,0:lx2+1,1:lx3))), &
+!              any(ieee_is_nan(J3halo(1:lx1,1:lx2,0:lx3+1))), &
+!              any(ieee_is_nan(divtmp(1:lx1,1:lx2,1:lx3)))
 !-------
 
 
@@ -623,4 +628,4 @@
 
 end procedure potential_root_mpi_curv
 
-end submodule potential_root
+end submodule potential_root

src/numerical/potential/potential_worker.f90

@@ -94,6 +94,10 @@
              J3halo(0:lx1+1,0:lx2+1,0:lx3+1),x,0,lx1+1,0,lx2+1,0,lx3+1)
 srcterm=divtmp(1:lx1,1:lx2,1:lx3)
 !-------
+!print*, myid, any(ieee_is_nan(J1halo(0:lx1+1,1:lx2,1:lx3))), &
+!              any(ieee_is_nan(J2halo(1:lx1,0:lx2+1,1:lx3))), &
+!              any(ieee_is_nan(J3halo(1:lx1,1:lx2,0:lx3+1))), &
+!              any(ieee_is_nan(divtmp(1:lx1,1:lx2,1:lx3)))
 
 
 !-------
@@ -480,4 +484,4 @@
 
 end procedure potential_workers_mpi
 
-end submodule potential_worker
+end submodule potential_worker