Testing altered ring vars, extra layer of xyz2ijk guessing

2026-01-09 15:17:56 -05:00 · 2020-01-14 13:59:59 -05:00
parent f51cb574a3
commit 3fa5f88b3c
8 changed files with 74 additions and 25 deletions
--- a/cmake/compilers.cmake
+++ b/cmake/compilers.cmake
@@ -51,7 +51,7 @@ if(CMAKE_Fortran_COMPILER_ID MATCHES Intel)
 	set(PROD "-align array64byte -align rec32byte -no-prec-div -fast-transcendentals")
 	#Debug
 	set(DEBUG "-g -traceback -check bounds -check uninit -debug all -gen-interfaces -warn interfaces -fp-stack-check")
-        set(PRODWITHDEBUGINFO "-O3 -g -traceback -debug all -align array64byte -align rec32byte -no-prec-div -fast-transcendentals")
+	set(PRODWITHDEBUGINFO "-O3 -g -traceback -debug all -align array64byte -align rec32byte -no-prec-div -fast-transcendentals")

 	#Now do OS-dep options
 	if (CMAKE_SYSTEM_NAME MATCHES Darwin)
@@ -59,15 +59,16 @@ if(CMAKE_Fortran_COMPILER_ID MATCHES Intel)
 	else()
 		#If we're not doing Mac, then add IPO
 		string(APPEND PROD " -ipo")
-                string(APPEND PRODWITHDEBUGINFO " -ipo")
+		string(APPEND PRODWITHDEBUGINFO " -ipo")
 	endif()

 	#Handle individual hosts
 	if (HOST MATCHES cheyenne)
 		string(APPEND PROD " -march=corei7 -axCORE-AVX2")
-                string(APPEND PRODWITHDEBUGINFO " -march=corei7 -axCORE-AVX2")
+		string(APPEND PRODWITHDEBUGINFO " -march=corei7 -axCORE-AVX2")
 	endif()
 	
+
 elseif(CMAKE_Fortran_COMPILER_ID MATCHES GNU)
 	set(dialect "-ffree-form -ffree-line-length-none -fimplicit-none")
 	#Base
--- a/src/chimp/gridloc.F90
+++ b/src/chimp/gridloc.F90
@@ -229,7 +229,7 @@ module gridloc

        logical :: isIn,inCell
        real(rp) :: xs,ys,lfmC(NDIM)
-        real(rp) :: xp(2), xCs(4,2)
+        real(rp) :: xp(2)
        integer :: i,j,i0,j0,k0
        integer :: i1,i2,j1,j2

@@ -253,26 +253,52 @@ module gridloc
        k0 = min(floor(lfmC(KDIR)/locAux%dTh) +1,ebGr%Nkp) !Evenly spaced k
        ijk(KDIR) = k0

+        xp = [xs,ys]
        !Use provided guess if present,
        if (present(ijkO)) then
-            xCs(1,:) = ebGr%xyz(ijkO(IDIR)  ,ijkO(JDIR)  ,ebGr%ks,XDIR:YDIR)
-            xCs(2,:) = ebGr%xyz(ijkO(IDIR)+1,ijkO(JDIR)  ,ebGr%ks,XDIR:YDIR)
-            xCs(3,:) = ebGr%xyz(ijkO(IDIR)+1,ijkO(JDIR)+1,ebGr%ks,XDIR:YDIR)
-            xCs(4,:) = ebGr%xyz(ijkO(IDIR),  ijkO(JDIR)+1,ebGr%ks,XDIR:YDIR)
-            xp = [xs,ys]
            !Test guess
-            inCell = inCell2D(xp,xCs)
-            if (inCell) then
-                !Found it, let's get out of here
-                ijk(IDIR) = ijkO(IDIR)
-                ijk(JDIR) = ijkO(JDIR)
+            inCell = CheckIJ(xp,ijkO(IDIR:JDIR),Model,ebGr)
+            if (inCell) then !Found it, let's get out of here
+                ijk(IDIR:JDIR) = ijkO(IDIR:JDIR)
+                !write(*,*) 'Found guess!'
                return
            endif
-        endif
+            !Wasn't in original guess, check adjacent cells
+            i1 = max(ijkO(IDIR)-1,ebGr%is)
+            i2 = min(ijkO(IDIR)+1,ebGr%ie)
+            j1 = max(ijkO(JDIR)-1,ebGr%js)
+            j2 = min(ijkO(JDIR)+1,ebGr%je)
+
+            do i=i1,i2
+                do j=j1,j2
+                    inCell = CheckIJ(xp,[i,j],Model,ebGr)
+                    if (inCell) then !Found it, let's get out of here
+                        ijk(IDIR:JDIR) = [i,j]
+                        !write(*,*) 'Found 1st halo!'
+                        return
+                    endif
+                enddo
+            enddo
+
+        endif !Using guess

        !If we're still here, do this the hard way
        !Cut out obviously incorrect 2D indices
        call lfmChop(Model,ebGr,[xs,ys],i1,i2,j1,j2)
+
+        ! !Start by looping
+        ! do i=i1,i2
+        !     do j=j1,j2
+        !         inCell = CheckIJ(xp,[i,j],Model,ebGr)
+        !         if (inCell) then !Found it, let's get out of here
+        !             ijk(IDIR:JDIR) = [i,j]
+        !             !write(*,*) 'Found in chop halo!'
+        !             return
+        !         endif
+        !     enddo
+        ! enddo
+
+        !If still here, just pick the closest one
        ijk(IDIR:JDIR) = minloc( (locAux%xxC(i1:i2,j1:j2)-xs)**2.0 + (locAux%yyC(i1:i2,j1:j2)-ys)**2.0 )
        ijk(IDIR:JDIR) = ijk(IDIR:JDIR) + [i1-1,j1-1] !Correct for offset
        if (ijk(KDIR)<0) then
@@ -281,6 +307,26 @@ module gridloc

    end subroutine Loc_LFM

+    !Check whether 2D point xy is in LFM cell ijG
+    function CheckIJ(xy,ijG,Model,ebGr) result(isIn)
+        real(rp), intent(in) :: xy (NDIM-1)
+        integer, intent(in)  :: ijG(NDIM-1)
+        type(chmpModel_T), intent(in) :: Model
+        type(ebGrid_T), intent(in) :: ebGr
+
+        logical :: isIn
+
+        real(rp) :: xCs(4,2)
+
+        xCs(1,:) = ebGr%xyz(ijG(IDIR)  ,ijG(JDIR)  ,ebGr%ks,XDIR:YDIR)
+        xCs(2,:) = ebGr%xyz(ijG(IDIR)+1,ijG(JDIR)  ,ebGr%ks,XDIR:YDIR)
+        xCs(3,:) = ebGr%xyz(ijG(IDIR)+1,ijG(JDIR)+1,ebGr%ks,XDIR:YDIR)
+        xCs(4,:) = ebGr%xyz(ijG(IDIR),  ijG(JDIR)+1,ebGr%ks,XDIR:YDIR)
+
+        !Test guess
+        isIn = inCell2D(xy,xCs)
+     
+    end function CheckIJ
    !3D localization routine for spherical grid
    subroutine Loc_SPH(xyz,ijk,Model,ebGr,isInO,ijkO)
        real(rp), intent(in) :: xyz(NDIM)
@@ -433,6 +479,7 @@ module gridloc
        j1 = max(j1-1,ebGr%js)
        j2 = min(j2+1,ebGr%je)

+        !write(*,*) 'Chop down to i1,i2,j1,j2 = ', i1,i2,j1,j2
    end subroutine lfmChop


--- a/src/drivers/voltronx.F90
+++ b/src/drivers/voltronx.F90
@@ -31,7 +31,7 @@ program voltronx
    !Do any updates to Voltron
        call stepVoltron(vApp,gApp)
        
-    !Coupling    
+    !Coupling
        call Tic("DeepCoupling")
        if ( (vApp%time >= vApp%DeepT) .and. vApp%doDeep ) then
            call DeepUpdate(vApp, gApp, vApp%time)
--- a/src/gamera/apps/wind.F90
+++ b/src/gamera/apps/wind.F90
@@ -320,8 +320,8 @@ module wind
        
        !$OMP PARALLEL DO default(shared) collapse(2) &
        !$OMP private(i,j,k,n,xcc,nHat,e1,e2,ecc,Vxyz,Bxyz,swExyz,wSW,D,P,mhdExyz)
-        do k=Grid%ksg,Grid%keg
-            do j=Grid%jsg,Grid%jeg
+        do k=Grid%ks,Grid%ke+1
+            do j=Grid%js,Grid%je+1
                do i=Grid%ie-2,Grid%ie+1
                    xcc = Grid%xyzcc(i,j,k,:)
                    nHat = Grid%Tf(ip+1,j,k,NORMX:NORMZ,IDIR)
--- a/src/gamera/gamapp.F90
+++ b/src/gamera/gamapp.F90
@@ -77,7 +77,7 @@ module gamapp
        call Tic("BCs")
        call EnforceBCs(gameraApp%Model,gameraApp%Grid,gameraApp%State)
        call Toc("BCs")
-
+        
    end subroutine stepGamera

    subroutine UpdateStateData(gameraApp)
--- a/src/gamera/mhdgroup.F90
+++ b/src/gamera/mhdgroup.F90
@@ -38,6 +38,7 @@ module mhdgroup
                call Toc("HackE")
            endif
            call Toc("E-Field")
+            
        endif
        
        !Get plasma stresses
--- a/src/gamera/ringutils.F90
+++ b/src/gamera/ringutils.F90
@@ -23,7 +23,7 @@ module ringutils
    !Which ring vars
    !doRAVarE = T: rho,mom,inte
    !doRAVarE = F: rho,mom,rho*Cs^2
-    logical, parameter :: doRAVarE = .true.
+    logical, parameter :: doRAVarE = .false.

    contains

--- a/src/gamera/stress.F90
+++ b/src/gamera/stress.F90
@@ -549,14 +549,14 @@ module stress
            By = MagLRB(i,YDIR)
            Bz = MagLRB(i,ZDIR)

+            B0x = B0(i,XDIR)
+            B0y = B0(i,YDIR)
+            B0z = B0(i,ZDIR)
+
            dPb = 0.5*(Bx**2.0 + By**2.0 + Bz**2.0) !Pressure in residual field
            if (Model%doBackground) then
-                B0x = B0(i,XDIR)
-                B0y = B0(i,YDIR)
-                B0z = B0(i,ZDIR)
-                Va2 = ( (Bx+B0x)**2.0 + (By+B0y)**2.0 + (Bz+B0z)**2.0)/D
+                Va2 = ( (Bx+B0x)**2.0 + (By+B0y)**2.0 + (Bz+B0z)**2.0 )/D
                bbD(i) = bbD(i) + 0.5*( (Bx+B0x)**2.0 + (By+B0y)**2.0 + (Bz+B0z)**2.0)
-
            else
                Va2 = 2*dPb/D
                bbD(i) = bbD(i) + 0.5*(Bx**2.0 + By**2.0 + Bz**2.0)