Skip to content

Commit

Permalink
Revert "Uses precomputed logtg"
Browse files Browse the repository at this point in the history 
This reverts commit 72fd093.
  • Loading branch information
@holm10
holm10 committed Jan 16, 2025
1 parent 72fd093 commit c8102df
Show file tree
Hide file tree
Showing 8 changed files with 319 additions and 212 deletions.
10 changes: 4 additions & 6 deletions bbb/bbb.v
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1287,8 +1287,6 @@ qe real /1.6022e-19/ #Elementary charge
mu0 real /1.25663706e-6/ #Vac. magnetic perm.
eps0 real /8.8542e-12/ #Vac. dielectric perm.
rt8opi real /1.595769121606e0/ #sqrt(8/pi)
log_10 real /2.302585092994046/ #log(10)
log10ev real /-18.795283272599516/ #log10(ev)

***** Comtra restart:
#Variables that contain the transport parameters.
Expand Down Expand Up @@ -1595,7 +1593,7 @@ mg(1:ngsp) _real [kg] /1.67e-27/ #gas species mass, calc. fr minu
facmg(1:nispmx) real /nispmx*1./ #scale factor for mg to recov old case
znucl(1:nisp) _integer [ ] #tot. nucl. charge, calc. from znuclin
ni(0:nx+1,0:ny+1,1:nisp) _real [m^-3] +threadprivate #ion density in primary cell (ix,iy)
logni(0:nx+1,0:ny+1,1:nisp) _real [m^-3] +threadprivate #log ion density in primary cell (ix,iy)
lni(0:nx+1,0:ny+1,1:nisp) _real [m^-3] #log(ion dens) in prim. cell (ix,iy)
nm(0:nx+1,0:ny+1,1:nisp) _real [kg*m^-3] +threadprivate #mass density [nm(,,1) is sum, exclud.
#gas, if nusp=1, isimpon=5] in cell
nz2(0:nx+1,0:ny+1) _real [m^-3] +threadprivate #sum of ni*zi**2 over all ion species
Expand Down Expand Up @@ -1643,6 +1641,7 @@ logte(0:nx+1,0:ny+1) _real [J] +threadprivate #log electron temperature
logti(0:nx+1,0:ny+1) _real [J] +threadprivate #log ion temperature in primary cell
ng(0:nx+1,0:ny+1,1:ngsp) _real [m^-3] +threadprivate #gas density in primary cell (ix,iy)
logng(0:nx+1,0:ny+1,1:ngsp) _real [m^-3] +threadprivate #log gas density in primary cell (ix,iy)
lng(0:nx+1,0:ny+1,1:ngsp) _real [m^-3] +threadprivate #log(gas dens) in prim. cell (ix,iy)
uug(0:nx+1,0:ny+1,1:ngsp) _real [m/s] +threadprivate #ratio gas-flux/density at x-face;
#if orthog mesh, poloidal gas velocity
uuxg(0:nx+1,0:ny+1,1:ngsp) _real [m/s] +threadprivate #poloidal-only component of uug;
Expand Down Expand Up @@ -3063,9 +3062,8 @@ ix2g(0:nx+1,0:ny+1) _integer #ix index for sec. interm. (ix,iy) pt.
iy2g(0:nx+1,0:ny+1) _integer #iy index for sec. interm. (ixo,iy) pt.
ixv2g(0:nx+1,0:ny+1) _integer #ixv index for sec. interm.(ixvo,iy) pt.
iyv2g(0:nx+1,0:ny+1) _integer #iyv index for sec.interm.(ixvo,iy) pt.
nis(0:nxold+1,0:nyold+1,1:nisp) _real [m^-3] +state
#ion dens at last success. calc
lognis(0:nxold+1,0:nyold+1,1:nisp) _real [m^-3] +restart #log ion dens at last success. calc
nis(0:nxold+1,0:nyold+1,1:nisp) _real [m^-3] +state
#ion dens at last success. calc
tes(0:nxold+1,0:nyold+1) _real [J] #elec. temp at last success. calc
tis(0:nxold+1,0:nyold+1) _real [J] #ion temp at last success. calc
tgs(0:nxold+1,0:nyold+1,1:ngsp) _real [J] #gas temp at last success. calc
Expand Down
4 changes: 2 additions & 2 deletions bbb/boundary.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4827,7 +4827,7 @@ c_mpi call mpi_allreduce(sumarr,gsumarr,2,MPI_DOUBLE_PRECISION,MPI_SUM,
c Compute fluxes along left (inner divertor - need to change sign)
do jx = 1, nxpt #number of x-points; =1 for single null
do iy = 1, ny
sxcos = sx(ixlb(jx),iy)*cosangfx(ixlb(jx),iy)
sxcos = sx(ixlb(jx),iy)*cos(angfx(ixlb(jx),iy))
ue_ion_flux_sum = 0.
do ifld = 1, 1 #placeholder; needs extension for impurities
ue_part_fluxh2p1lb(iy,jx) = -fnix(ixlb(jx),iy,1)/sxcos
Expand Down Expand Up @@ -4858,7 +4858,7 @@ c Compute fluxes along right (outer divertor)
nxt = ixrb(jx)-1 #effec nx
nxt1 = ixrb(jx) #effec nx+1
do iy = 1, ny
sxcos = sx(nxt,iy)*cosangfx(nxt,iy)
sxcos = sx(nxt,iy)*cos(angfx(nxt,iy))
ue_ion_flux_sum = 0.
do ifld = 1, 1 #placeholder; needs extension for impurities
ue_part_fluxh2p1rb(iy,jx) = fnix(nxt,iy,1)/sxcos
Expand Down
20 changes: 6 additions & 14 deletions bbb/convert.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@
if(ineudif .ne. 3) then
yl(iv) = ng(ix,iy,igsp)/n0g(igsp)
elseif(ineudif .eq. 3) then
yl(iv) = logng(ix,iy,igsp)
yl(iv) = lng(ix,iy,igsp)
endif
rtol(iv) = rtolv(igrid)*bfac
atol(iv) = cngatol*rtol(iv)*bfac*abs(yl(iv))
Expand Down Expand Up @@ -261,8 +261,7 @@ c_mpi integer status(MPI_STATUS_SIZE)
ne(ix,iy) = ne(ix,iy) + zi(ifld)*ni(ix,iy,ifld)
if (isupgon(1).eq.1 .and. zi(ifld).eq.0) then
ng(ix,iy,1) = ni(ix,iy,ifld)
logng(ix,iy,1)=log(abs(ng(ix,iy,1)))
if (ineudif .eq. 3) logng(ix,iy,1)=log(ng(ix,iy,1))
if (ineudif .eq. 3) lng(ix,iy,1)=log(ng(ix,iy,1))
else
nit(ix,iy) = nit(ix,iy) + ni(ix,iy,ifld)
if (isimpon.ge.5 .and. nusp_imp.eq.0)
Expand All @@ -281,7 +280,6 @@ c_mpi integer status(MPI_STATUS_SIZE)
if(isteonxy(ix,iy) .eq. 1) then
te(ix,iy)=yl(idxte(ix,iy))*ennorm/(1.5*ntemp)
te(ix,iy) = max(te(ix,iy), temin*ev) #NEW Feb4,2018
logte(ix,iy)=log(abs(te(ix,iy)))
endif
do 65 igsp =1, ngsp
if(isngonxy(ix,iy,igsp) .eq. 1) then
Expand All @@ -294,10 +292,9 @@ c_mpi integer status(MPI_STATUS_SIZE)
igspneg = igsp
endif
elseif(ineudif .eq. 3) then
logng(ix,iy,igsp) = yl(idxg(ix,iy,igsp))
ng(ix,iy,igsp) = exp(logng(ix,iy,igsp))
lng(ix,iy,igsp) = yl(idxg(ix,iy,igsp))
ng(ix,iy,igsp) = exp(lng(ix,iy,igsp))
endif
logng(ix,iy,igsp)=log(abs(ng(ix,iy,igsp)))
endif
if(istgonxy(ix,iy,igsp) .eq. 1) then
ntemp = ng(ix,iy,igsp)
Expand All @@ -306,7 +303,6 @@ c_mpi integer status(MPI_STATUS_SIZE)
. (1.5*ntemp)
tg(ix,iy,igsp) = max(tg(ix,iy,igsp), tgmin*ev)
endif
logtg(ix,iy,igsp)=log(abs(tg(ix,iy,igsp)))
65 continue
ntemp = nit(ix,iy) + cngtgx(1)*ng(ix,iy,1)
if(isflxvar .eq. 0) ntemp = nnorm
Expand Down Expand Up @@ -519,11 +515,8 @@ subroutine convsr_aux (ixl, iyl)
do 12 iy = js, je
do 11 ix = is, ie
pri(ix,iy,ifld) = ni(ix,iy,ifld) * ti(ix,iy)
logpri(ix,iy,ifld) = logni(ix,iy,ifld) + logti(ix,iy)
if (ifld.eq.iigsp .and. istgon(1).eq.1) then
pri(ix,iy,ifld) = ni(ix,iy,ifld) * tg(ix,iy,1)
logpri(ix,iy,ifld) = logni(ix,iy,ifld) + logtg(ix,iy,1)
endif
if (ifld.eq.iigsp .and. istgon(1).eq.1) pri(ix,iy,ifld) =
. ni(ix,iy,ifld) * tg(ix,iy,1)
if (zi(ifld).ne.0.) then
pr(ix,iy) = pr(ix,iy) + pri(ix,iy,ifld)
zeff(ix,iy)=zeff(ix,iy)+zi(ifld)**2*ni(ix,iy,ifld)
Expand Down Expand Up @@ -559,7 +552,6 @@ ccc ne(ix,iy) = rnec (ni(ix,iy,1), nzsp, ni(ix,iy,nhsp+1),
. (1-istgcon(igsp))*rtg2ti(igsp)*ti(ix,iy) +
. istgcon(igsp)*tgas(igsp)*ev
pg(ix,iy,igsp) = ng(ix,iy,igsp)*tg(ix,iy,igsp)
logpg(ix,iy,igsp) = logng(ix,iy,igsp) + logtg(ix,iy,igsp)
enddo
15 continue
16 continue
Expand Down
34 changes: 5 additions & 29 deletions bbb/geometry.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -868,13 +868,13 @@ call s2fill (nx+2, ny+2, 0., fxpyv(0:nx+1,0:ny+1,iu), 1, nx+2)
call s2fill (nx+2, ny+2, 0., fymxv(0:nx+1,0:ny+1,iu), 1, nx+2)
call s2fill (nx+2, ny+2, 0., fypxv(0:nx+1,0:ny+1,iu), 1, nx+2)
enddo
if (isnonog .ge. 1) then
call nonorthg # Sets dist btwn interp pts as dxnog, dynog
if (isnonog .ge. 1) then
call nonorthg # Sets dist btwn interp pts as dxnog, dynog
# Also sets nonorth stencils:fx0,fxm,fy0,fym etc
do ix = 0, nx+1 # Bdry value set: matters little except vygtan
ccc gxfn(ix,0) = gxfn(ix,1)
dxnog(ix,0) = dxnog(ix,1)
dxnog(ix,ny+1) = dxnog(ix,ny)
dxnog(ix,ny+1) = dxnog(ix,ny)
enddo
do iy = 0, ny+1 # likewise, dxnog(nx+1,) not relevant, but avoid 0
dxnog(nx+1,iy) = 0.1*dxnog(nx,iy)
Expand Down Expand Up @@ -906,10 +906,10 @@ ccc gxfn(ix,0) = gxfn(ix,1)
dx(ix,iy) = max(dx(ix,iy), dxmin)
cc for nonorthogonal mesh, dy gets reduced by cosine(average_angle)
if ((islimon .ne. 0) .and. (ix .eq. ix_lim+1)) then
dy(ix,iy) = dy(ix,iy)*cos(angfx(ix,iy))
dy(ix,iy) = dy(ix,iy)*cos( angfx(ix,iy) )
elseif (nxpt==2 .and. ix==ixlb(nxpt)) then # full double null
c Effectively, use angfx(ix-1,iy)=angfx(ix,iy) at left boundaries
dy(ix,iy) = dy(ix,iy)*cos(angfx(ix,iy))
dy(ix,iy) = dy(ix,iy)*cos( angfx(ix,iy) )
else
dy(ix,iy) = dy(ix,iy)*cos( 0.5*
. (angfx(ix,iy) + angfx(ixm1(ix,iy),iy)) )
Expand Down Expand Up @@ -1957,11 +1957,9 @@ c of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
iym1 = max(0, iy-1)
do ix = 0, nx+1
angfx(ix,iy) = 0.5*(vtag(ix,iy)+vtag(ix,iym1))
cosangfx(ix,iy) = COS(angfx(ix,iy))
if (islimon .eq. 1) then
if (iy.eq.iy_lims) then
angfx(ix,iy) = vtag(ix,iy)
cosangfx(ix,iy) = COS(angfx(ix,iy))
endif
if (iy .eq. iy_lims-1) then
angfx(ix,iy) = vtag(ix,iym1)
Expand All @@ -1978,49 +1976,38 @@ c of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
if(redopltvtag == 1) then
do iy = 0, ny+1
angfx(0,iy) = 2*angfx(1,iy) - angfx(2,iy)
cosangfx(0,iy) = COS(angfx(0,iy))
vtag(0,iy) = 2*vtag(1,iy) - vtag(2,iy)
if (nxomit .gt. 0) then
angfx(0,iy) = angfx(1,iy)
vtag(0,iy) = vtag(1,iy)
endif
if (islimon .ne. 0) then
angfx(ix_lim-1,iy) = 2*angfx(ix_lim-2,iy)-angfx(ix_lim-3,iy)
cosangfx(ix_lim-1,iy) = COS(angfx(ix_lim-1,iy))
vtag(ix_lim-1,iy) = 2*vtag(ix_lim-2,iy)-vtag(ix_lim-3,iy)
angfx(ix_lim,iy) = angfx(ix_lim-1,iy)
cosangfx(ix_lim,iy) = COS(angfx(ix_lim,iy))
vtag(ix_lim,iy) = vtag(ix_lim-1,iy)
angfx(ix_lim+1,iy) = 2*angfx(ix_lim+2,iy)-angfx(ix_lim+3,iy)
cosangfx(ix_lim+1,iy) = COS(angfx(ix_lim+1,iy))
vtag(ix_lim+1,iy) = 2*vtag(ix_lim+2,iy)-vtag(ix_lim+3,iy)
endif
if (nxpt==2) then # full double null
angfx(ixrb(1),iy) = 2*angfx(ixrb(1)-1,iy) - angfx(ixrb(1)-2,iy)
cosangfx(ixrb(1),iy) = COS(angfx(ixrb(1),iy))
vtag(ixrb(1),iy) = 2*vtag(ixrb(1)-1,iy) - vtag(ixrb(1)-2,iy)
angfx(ixrb(1)+1,iy) = angfx(ixrb(1),iy)
cosangfx(ixrb(1)+1,iy) = COS(angfx(ixrb(1)+1,iy))
vtag(ixrb(1)+1,iy) = vtag(ixrb(1),iy)
angfx(ixlb(2),iy) = 2*angfx(ixlb(2)+1,iy) - angfx(ixlb(2)+2,iy)
cosangfx(ixlb(2),iy) = COS(angfx(ixlb(2),iy))
vtag(ixlb(2),iy) = 2*vtag(ixlb(2)+1,iy) - vtag(ixlb(2)+2,iy)
endif
angfx(nx,iy) = 2*angfx(nx-1,iy) - angfx(nx-2,iy)
cosangfx(nx,iy) = COS(angfx(nx,iy))
vtag(nx,iy) = 2*vtag(nx-1,iy) - vtag(nx-2,iy)
angfx(nx+1,iy) = angfx(nx,iy) # not really used
cosangfx(nx+1,iy) = COS(angfx(nx+1,iy))
vtag(nx+1,iy) = vtag(nx,iy) # not really used
enddo
endif

c... Set angfx and vtag in the y-guard cells to adjacent values
do ix = 0, nx+1
angfx(ix,0) = angfx(ix,1)
cosangfx(ix,0) = COS(angfx(ix,0))
angfx(ix,ny+1) = angfx(ix,ny)
cosangfx(ix,ny+1) = COS(angfx(ix,ny+1))
vtag(ix,0) = vtag(ix,1)
vtag(ix,ny+1) = vtag(ix,ny)
enddo
Expand All @@ -2030,33 +2017,25 @@ c of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
do jx = 1, nxpt
if (ixpt1(jx) .gt. 0) then
angfx(ixpt1(jx),iysptrx1(jx)-1) = 0.
cosangfx(ixpt1(jx),iysptrx1(jx)-1) = COS(0.)
angfx(ixpt1(jx),iysptrx1(jx) ) = 0.
cosangfx(ixpt1(jx),iysptrx1(jx)) = COS(0.)
angfx(ixpt1(jx),iysptrx1(jx)+1) = 0.
cosangfx(ixpt1(jx),iysptrx1(jx)+1) = COS(0.)
endif
if (ixpt2(jx) .gt. 0) then
angfx(ixpt2(jx),iysptrx2(jx)-1) = 0.
cosangfx(ixpt2(jx),iysptrx2(jx)-1) = COS(0.)
angfx(ixpt2(jx),iysptrx2(jx) ) = 0.
cosangfx(ixpt2(jx),iysptrx2(jx)) = COS(0.)
angfx(ixpt2(jx),iysptrx2(jx)+1) = 0.
cosangfx(ixpt2(jx),iysptrx2(jx)+1) = COS(0.)
endif
enddo

c... reset values next to cut at ixpt1 or ixpt2 if half-space problem
if (isfixlb(1).eq.2 .and. ixpt2(1).gt.0) then
do iy = 0, iysptrx1(1)
angfx(ixpt2(1),iy) = 0.
cosangfx(ixpt2(1),iy) = COS(0.)
enddo
endif
if (isfixrb(1).eq.2 .and. ixpt1(1).gt.0) then
do iy = 0, iysptrx1(1)
angfx(ixpt1(1),iy) = 0.
cosangfx(ixpt1(1),iy) = COS(0.)
enddo
endif

Expand All @@ -2065,11 +2044,8 @@ c of cell (ix,iy) and vertex 1 of cell (ix+1,iy+1).
if ((isudsym==1.or.(geometry.eq.'dnXtarget')) .and. nxc.gt.0) then
do iy = 0, ny+1
angfx(nxc-1,iy) = 0.
cosangfx(nxc-1,iy) = COS(0.)
angfx(nxc ,iy) = 0.
cosangfx(nxc,iy) = COS(0.)
angfx(nxc+1,iy) = 0.
cosangfx(nxc+1,iy) = COS(0.)
vtag(nxc-1,iy) = 0.
vtag(nxc ,iy) = 0.
vtag(nxc+1,iy) = 0.
Expand Down
Loading

0 comments on commit c8102df

Please sign in to comment.