Bugfix/snowflake gt90

bbb/boundary.m

@@ -106,14 +106,23 @@ c_mpi      Use(MpiVars)  #module defined in com/mpivarsmod.F.in
 c...  do the iy = 0 boundary
 c...  if extrapolation b.c. on p.f. region, isextrpf=1, otherwise isextrpf=0
       if (iymnbcl .eq. 0) goto 1100   # interior domain boundary; no bdry eqn
-      ixc1 = max(0, ixpt1(1)+1)       # 1st  core cell;used for core flux BC
-      if ((geometry(1:9)=="snowflake" .and. geometry(10:11).ne."15")
+
+      if (geometry=='snowflake135' .or. geometry=='snowflake105') then 
+        ixc1 = max(0, ixpt2(1)+1) # 1st  core cell;used for core flux BC
+      else 
+        ixc1 = max(0, ixpt1(1)+1) # 1st  core cell;used for core flux BC
+      end if
+      if (geometry=='snowflake135' .or. geometry=='snowflake105') then 
+        ix_fl_bc = min(ixpt1(2), nx)
+      else if ((geometry(1:9)=="snowflake" .and. geometry(10:11).ne."15")
      &                     .or. geometry=="dnXtarget") then
-        ix_fl_bc = min(ixpt2(1), nx)
+        ix_fl_bc = min(ixpt2(1), nx) # last core cell;use for flux BC
       else
         ix_fl_bc = min(ixpt2(nxpt), nx) # last core cell;use for flux BC
       endif 
 
+      ! write(*,*) ix_fl_bc
+
 c  Note: j3 is local range index for iy passed from pandf in oderhs.m
       if (j3 .le. isextrnpf .or. j3 .le. isextrtpf .or. 
      .                                         j3 .le. isextrngc) then
@@ -227,7 +236,7 @@ cc                  if (fng_chem .ne. 0.) yldot(iv1) = -nurlxg*(
                      yldot(iv1) = -nurlxn*( ni(ix,0,ifld) -
      .                                ni(ixp1(ix,0),0,ifld) ) / n0(ifld)
                      if (ix.eq.ix_fl_bc) then
-                        ii=max(0,ixpt1(1)+1)
+                        ii=ixc1
                         fniytotc=fniy(ii,0,ifld)-fniycbo(ii,ifld)
                         fngytotc = fngy(ii,0,1)  # needs generalization
                         do
@@ -333,7 +342,7 @@ call xerrab ('** isnicore value not valid option **')
      .                          (rm(ixp1(ix,0),0,2)*ni(ix,0,ifld)) )/vpnorm
                     if (ix.eq.ix_fl_bc) then # not all Jac elems included
                       # sum core parallel momentum flux:
-                      ii = max(0, ixpt1(1)+1)
+                      ii = ixc1
                       fmiytotc = rm(ii,0,0)*fmiy(ii,0,ifld)
                       sytotc = sy(ii,0)
                       uztotc = uz(ii,0,ifld)/gxf(ii,0)
@@ -413,7 +422,7 @@ c     Force corner cells (0,0) and (nx+1,0) to relax to an average of
                      yldot(iv1) = nurlxn*( ni(ix,0,iimp) -
      .                                 ni(ixp1(ix,0),0,iimp) ) / n0(iimp)
                      if (ix.eq.ix_fl_bc) then
-                        ii=max(0,ixpt1(1)+1)
+                        ii=ixc1
                         fniytotc=fniy(ii,0,iimp)-fniycbo(ii,iimp)
                         do
                           ii=ixp1(ii,0)
@@ -488,7 +497,7 @@ c     Force corner cells (0,0) and (nx+1,0) to relax to an average of
         feiytotc = 0.
         do ix = i4+1-ixmnbcl, i8-1+ixmxbcl
           if (iymnbcl==1 .and. isixcore(ix)==1) then  # domain part of core bdry
-            ii = max(0, ixpt1(1)+1)
+            ii = ixc1
             feeytotc = feey(ii,0)-feeycbo(ii)
             feiytotc = feiy(ii,0)-feiycbo(ii)
             do    # loop over ii as changed by ii=ixp1 statement
@@ -530,7 +539,7 @@ c_mpi      call mpi_allreduce(feiytotc, sfeiytotc, 1,
                  yldot(iv1) = -nurlxe*(te(ix,0)-te(ixp1(ix,0),0))*n0(1)/ennorm
                  if (ix.eq.ix_fl_bc) then # not all Jac elems included
                     # sum core power:
-                    ii = max(0, ixpt1(1)+1)
+                    ii = ixc1
                     feeytotc = feey(ii,0)-feeycbo(ii)
                     do    # loop over ii as changed by ii=ixp1 statement
                        ii = ixp1(ii,0)
@@ -578,7 +587,7 @@ c_mpi      call mpi_allreduce(feiytotc, sfeiytotc, 1,
      .                                                     n0(1) /ennorm
                  if (ix.eq.ix_fl_bc) then # not all Jac elems included
                     # sum core power:
-                    ii = max(0, ixpt1(1)+1)
+                    ii = ixc1
                     feiytotc = feiy(ii,0)-feiycbo(ii)
                     do
                        ii = ixp1(ii,0)
@@ -990,7 +999,7 @@ call xerrab("***Input error: invalid istgpfc ***")
 c...  If isnewpot=1, phi boundary involves two eqns at iy=0 and 1
 c...  Note: j3 is local range index for iy passed from pandf in oderhs.m
       if (isnewpot*isphion.eq.1 .and. j3.le.3) then
-        do ix = min(i4+1-ixmnbcl,ixpt1(1)+1), max(i8-1+ixmxbcl,ixpt2(nxpt))
+        do ix = min(i4+1-ixmnbcl,ixc1), max(i8-1+ixmxbcl,ix_fl_bc)
           if(isphionxy(ix,0)*isphionxy(ix,1)==1) then
             iv  = idxphi(ix,0)
             iv1 = idxphi(ix,1)

bbb/geometry.m

@@ -1114,7 +1114,7 @@ ccc              sx(ix,iy) = xgbx*sx(ix,iy)
 
 c...  Calculate a normalization constant for the iy=0 cells
 c...  of the core boundary region
-      if ((geometry(1:9)=="snowflake" .and. geometry(10:11).ne."15")
+      if ((geometry(1:9)=="snowflake" .and. ((geometry(10:11).ne."15") .and. (geometry(10:12).ne."165")))
      &       .or. geometry=="dnXtarget") then
         ix_last_core_cell = ixpt2(1)
       else
@@ -1124,14 +1124,21 @@ ccc              sx(ix,iy) = xgbx*sx(ix,iy)
       ix = min(ix, nx+ixmxbcl)
       sygytotc = sy(ix,0)*gyf(ix,0)
       do ix = ixpt1(1)+1, ix_last_core_cell
+         !TODO: Should this block move below definition of isixcore below, and then this loop should only be over cells where isixcore(ix) == 1?
          if ( .not. (isudsym==1 .and. ((ix==nxc) .or. (ix==nxc+1))) ) then
-           sygytotc = sygytotc + sy(ix,0)*gyf(ix,0)
+         sygytotc = sygytotc + sy(ix,0)*gyf(ix,0)
          endif
       enddo
 
 c...  Setup the isixcore(ix) array: =1 if ix on iy=0 core bdry; =0 if not
       do ix = 0, nx+1
-        if ((geometry(1:9)=="snowflake" .and. geometry(10:11).ne."15")
+        if (geometry(1:9)=="snowflake" .and. (geometry(10:12)=="135" .or. geometry(10:12)=="105")) then 
+          if( ix > ixpt2(1) .and. ix <= ixpt1(2)) then
+            isixcore(ix) = 1
+          else
+            isixcore(ix) = 0
+          endif
+        else if ((geometry(1:9)=="snowflake" .and. ((geometry(10:11).ne."15") .and. (geometry(10:12).ne."165")))
      &       .or. geometry=="dnXtarget") then
           if( ix > ixpt1(1) .and. ix <= ixpt2(1)) then
             isixcore(ix) = 1
@@ -1149,7 +1156,6 @@ ccc              sx(ix,iy) = xgbx*sx(ix,iy)
           endif
         endif
       enddo
-
 *----------------------------------------------------------------------
 *  -- Calculate geometrical factors needed for curvature and grad_B drifts
 
@@ -1322,18 +1328,39 @@ ccc            gradby(ix,iy) = 0.
 *  -- define y on density faces -- at outboard midplane (?)
       if (ixpt2(1) > 0 .and. (isudsym.ne.1) .and.isddcon==0) then
          rmmax = rm(nxleg(1,1)+nxcore(1,1)+1,ny,0)
-         do ix = nxleg(1,1)+nxcore(1,1)+1, ixpt2(1)
-           if (rm(ix,ny,0) >= rmmax) then
-              rmmax = rm(ix,ny,0)
-              ixmp = ix
-           endif
-         enddo
+         if (geometry=='snowflake105' .or. geometry=='snowflake135') then 
+            do ix = max(0, ixpt2(1)+1), min(ixpt1(2), nx)
+               if (rm(ix,ny,0) >= rmmax) then
+                  rmmax = rm(ix,ny,0)
+                  ixmp = ix
+               endif
+            enddo
+         else if (geometry=='snowflake15') then 
+            do ix = max(0, ixpt1(1)+1) , min(ixpt2(2), nx)
+               if (rm(ix,ny,0) >= rmmax) then
+                  rmmax = rm(ix,ny,0)
+                  ixmp = ix
+               endif
+            enddo
+         else if (geometry(1:9)=="snowflake") then 
+            do ix = max(0, ixpt1(1)+1) , min(ixpt2(1), nx)
+               if (rm(ix,ny,0) >= rmmax) then
+                  rmmax = rm(ix,ny,0)
+                  ixmp = ix
+               endif
+            enddo
+         else 
+            do ix = nxleg(1,1)+nxcore(1,1)+1, ixpt2(1)
+            if (rm(ix,ny,0) >= rmmax) then
+               rmmax = rm(ix,ny,0)
+               ixmp = ix
+            endif
+            enddo
+         end if
       endif
       if (geometry.eq.'dnbot') ixmp = nxc+1
-      if (geometry.eq.'dnull' .or. geometry=='snowflake15' .or.
-     .    geometry=='snowflake45' .or. geometry=='snowflake75' .or.
+      if (geometry.eq.'dnull' .or.
      .    geometry=='dnXtarget' .or. geometry=='isoleg') ixmp = ixmdp(2)
-
 *  -- redefine ixmp if it is outside ix=0,nx domain
       if(ixmp.lt.0 .or. ixmp.gt.nx) then  #search for max rm
 	 rmmax = rm(nxomit,0,0)
@@ -1344,7 +1371,6 @@ ccc            gradby(ix,iy) = 0.
             endif
          enddo
       endif
-
 c     MER NOTE:
 c     The general case has multiple separatrices so yyf=0 is not
 c     uniquely defined; below we choose the innermost separatrix.

grd/grdread.m

@@ -170,6 +170,8 @@ call xerrab("**** requested grid data file not found")
 
       if (geometry=="dnull" .or. geometry=="snowflake15" .or.
      .    geometry=="snowflake45" .or. geometry=="snowflake75" .or.
+     .    geometry=="snowflake105" .or. geometry=="snowflake135" .or.
+     .    geometry=="snowflake165" .or.
      .    geometry=="dnXtarget" .or. geometry=="isoleg") then
          read(nuno,1999) nxm,nym
          read(nuno,1999) iysptrx1(1),iysptrx2(1)
@@ -224,6 +226,8 @@ call xerrab("**** requested grid data file not found")
 
       if (geometry=="dnull" .or. geometry=="snowflake15" .or.
      .    geometry=="snowflake45" .or. geometry=="snowflake75" .or.
+     .    geometry=="snowflake105" .or. geometry=="snowflake135" .or.
+     .    geometry=="snowflake165" .or.	  
      .    geometry=="dnXtarget" .or. geometry=="isoleg") then
          read(nuno,1999) nxm,nym
          read(nuno,1999) iysptrx1(1),iysptrx2(1)

src/uedge/gridue.py

@@ -6,8 +6,9 @@
 
 def read_gridpars(fname=None):
     from uedge import bbb, com, grd
-    dblxpts = ['dnull', 'snowflake15', 'snowflake45', 'snowflake75', 
-        'dnXtarget', 'isoleg']
+    dblxpts =   ['dnull', 'dnXtarget', 'isoleg'] + \
+                [f'snowflake{15+i*30}' for i in range(6)]
+
     geometry = com.geometry[0].decode('UTF-8').strip()
     if fname is None:
         fname = bbb.GridFileName[0].decode('UTF-8').strip()