diff --git a/bbb/boundary.m b/bbb/boundary.m index 094ab8a0..cf7216e2 100755 --- a/bbb/boundary.m +++ b/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) diff --git a/bbb/geometry.m b/bbb/geometry.m index e136b025..a10f97e1 100755 --- a/bbb/geometry.m +++ b/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. diff --git a/grd/grdread.m b/grd/grdread.m index c5944a76..dcffa8a4 100755 --- a/grd/grdread.m +++ b/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) diff --git a/src/uedge/gridue.py b/src/uedge/gridue.py index a823f152..831e663f 100644 --- a/src/uedge/gridue.py +++ b/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()