From e5942452b2f18b0aaa2d9d9d543f9255874a088c Mon Sep 17 00:00:00 2001 From: Matthew Zettergren Date: Tue, 12 May 2020 14:56:29 -0400 Subject: [PATCH] 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 | 6 +++--- src/numerical/calculus/gradient.f90 | 4 ++-- src/numerical/potential/potential_comm_mumps.f90 | 3 ++- src/numerical/potential/potential_root.f90 | 7 ++++++- src/numerical/potential/potential_worker.f90 | 6 +++++- 5 files changed, 18 insertions(+), 8 deletions(-) diff --git a/src/numerical/calculus/div.f90 b/src/numerical/calculus/div.f90 index 73caee687..7fbde26a0 100644 --- a/src/numerical/calculus/div.f90 +++ b/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 \ No newline at end of file +end submodule div diff --git a/src/numerical/calculus/gradient.f90 b/src/numerical/calculus/gradient.f90 index 2eb1619ae..c47570240 100644 --- a/src/numerical/calculus/gradient.f90 +++ b/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 diff --git a/src/numerical/potential/potential_comm_mumps.f90 b/src/numerical/potential/potential_comm_mumps.f90 index 7ec518916..a8c48e790 100644 --- a/src/numerical/potential/potential_comm_mumps.f90 +++ b/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 diff --git a/src/numerical/potential/potential_root.f90 b/src/numerical/potential/potential_root.f90 index f4094cf08..c89919e8f 100644 --- a/src/numerical/potential/potential_root.f90 +++ b/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 \ No newline at end of file +end submodule potential_root diff --git a/src/numerical/potential/potential_worker.f90 b/src/numerical/potential/potential_worker.f90 index 77bd52b6b..39404494b 100644 --- a/src/numerical/potential/potential_worker.f90 +++ b/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 \ No newline at end of file +end submodule potential_worker