Skip to content

Commit

Permalink
Uses precomputed logtg
Browse files Browse the repository at this point in the history 
Pre-computes COS(angfx)

- Removes obsolete subroutine coneq
- Sets up precomputed log variables

Removes logangfx usage in bbb.geometry.m

As it seems to have been causing issues with defining com.sx: not a big deal since geometry only needs to be calculated once

Uses precomputed logng

- Replaces lng
  • Loading branch information
@holm10
holm10 committed Jan 16, 2025
1 parent 7ae9248 commit 72fd093
Show file tree
Hide file tree
Showing 8 changed files with 212 additions and 319 deletions.
10 changes: 6 additions & 4 deletions bbb/bbb.v
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1287,6 +1287,8 @@ 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 @@ -1593,7 +1595,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)
lni(0:nx+1,0:ny+1,1:nisp) _real [m^-3] #log(ion dens) in prim. cell (ix,iy)
logni(0:nx+1,0:ny+1,1:nisp) _real [m^-3] +threadprivate #log ion density in primary 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 @@ -1641,7 +1643,6 @@ 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 @@ -3062,8 +3063,9 @@ 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
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
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)*cos(angfx(ixlb(jx),iy))
sxcos = sx(ixlb(jx),iy)*cosangfx(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)*cos(angfx(nxt,iy))
sxcos = sx(nxt,iy)*cosangfx(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: 14 additions & 6 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) = lng(ix,iy,igsp)
yl(iv) = logng(ix,iy,igsp)
endif
rtol(iv) = rtolv(igrid)*bfac
atol(iv) = cngatol*rtol(iv)*bfac*abs(yl(iv))
Expand Down Expand Up @@ -261,7 +261,8 @@ 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)
if (ineudif .eq. 3) lng(ix,iy,1)=log(ng(ix,iy,1))
logng(ix,iy,1)=log(abs(ng(ix,iy,1)))
if (ineudif .eq. 3) logng(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 @@ -280,6 +281,7 @@ 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 @@ -292,9 +294,10 @@ c_mpi integer status(MPI_STATUS_SIZE)
igspneg = igsp
endif
elseif(ineudif .eq. 3) then
lng(ix,iy,igsp) = yl(idxg(ix,iy,igsp))
ng(ix,iy,igsp) = exp(lng(ix,iy,igsp))
logng(ix,iy,igsp) = yl(idxg(ix,iy,igsp))
ng(ix,iy,igsp) = exp(logng(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 @@ -303,6 +306,7 @@ 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 @@ -515,8 +519,11 @@ 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)
if (ifld.eq.iigsp .and. istgon(1).eq.1) pri(ix,iy,ifld) =
. ni(ix,iy,ifld) * tg(ix,iy,1)
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 (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 @@ -552,6 +559,7 @@ 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: 29 additions & 5 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,9 +1957,11 @@ 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 @@ -1976,38 +1978,49 @@ 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 @@ -2017,25 +2030,33 @@ 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 @@ -2044,8 +2065,11 @@ 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 72fd093

Please sign in to comment.