Skip to content

Commit

Permalink
Uses precomputed logng
Browse files Browse the repository at this point in the history 
- Replaces lng
  • Loading branch information
@holm10
holm10 committed Jan 15, 2025
1 parent d242e90 commit ee50085
Show file tree
Hide file tree
Showing 5 changed files with 57 additions and 51 deletions.
1 change: 0 additions & 1 deletion bbb/bbb.v
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1643,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
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
80 changes: 40 additions & 40 deletions bbb/oderhs.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5832,16 +5832,16 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
. (ix==ixrb(jx).and.ixmxbcl==1) ) isxyfl = .false.
enddo
if (methgx .eq. 6) then # log interpolation
grdnv =( ( fym (ix,iy,1)*log(ng(ix2,iy1 ,igsp)) +
. fy0 (ix,iy,1)*log(ng(ix2,iy ,igsp)) +
. fyp (ix,iy,1)*log(ng(ix2,iy+1,igsp)) +
. fymx(ix,iy,1)*log(ng(ix ,iy1 ,igsp)) +
. fypx(ix,iy,1)*log(ng(ix, iy+1,igsp)) )
. -( fym (ix,iy,0)*log(ng(ix ,iy1 ,igsp)) +
. fy0 (ix,iy,0)*log(ng(ix ,iy ,igsp)) +
. fyp (ix,iy,0)*log(ng(ix ,iy+1,igsp)) +
. fymx(ix,iy,0)*log(ng(ix4,iy1 ,igsp)) +
. fypx(ix,iy,0)*log(ng(ix6,iy+1,igsp)) ) )/
grdnv =( ( fym (ix,iy,1)*logng(ix2,iy1 ,igsp) +
. fy0 (ix,iy,1)*logng(ix2,iy ,igsp) +
. fyp (ix,iy,1)*logng(ix2,iy+1,igsp) +
. fymx(ix,iy,1)*logng(ix ,iy1 ,igsp) +
. fypx(ix,iy,1)*logng(ix, iy+1,igsp) )
. -( fym (ix,iy,0)*logng(ix ,iy1 ,igsp) +
. fy0 (ix,iy,0)*logng(ix ,iy ,igsp) +
. fyp (ix,iy,0)*logng(ix ,iy+1,igsp) +
. fymx(ix,iy,0)*logng(ix4,iy1 ,igsp) +
. fypx(ix,iy,0)*logng(ix6,iy+1,igsp) ) )/
. dxnog(ix,iy)
elseif (methgx .eq. 7) then # inverse interpolation
grdnv =( 1/(fym (ix,iy,1)/ng(ix2,iy1 ,igsp) +
Expand Down Expand Up @@ -5874,9 +5874,9 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
. 0.5*(nuiz(ix,iy,igsp)+nuiz(ix2,iy,igsp))
if (methgx .eq. 6) then
fngxy(ix,iy,igsp) = exp( 0.5*
. (log(ng(ix2,iy,igsp))+log(ng(ix,iy,igsp))) )*
. (logng(ix2,iy,igsp)+logng(ix,iy,igsp)) )*
. difgx2*(grdnv/cosangfx(ix,iy) -
. (log(ng(ix2,iy,igsp)) - log(ng(ix,iy,igsp)))*
. (logng(ix2,iy,igsp) - logng(ix,iy,igsp))*
. gxf(ix,iy) ) * sx(ix,iy)
else
fngxy(ix,iy,igsp) = difgx2*( grdnv/cosangfx(ix,iy) -
Expand Down Expand Up @@ -6766,7 +6766,7 @@ c physically meaningful gas variables (flux, velocity, etc)
if (isupgon(igsp) .eq. 1) then
qtgf = qtgf + rrv(ix,iy)*up(ix,iy,iigsp)*sx(ix,iy)
endif
qsh = csh * (lng(ix,iy,igsp)-lng(ix2,iy,igsp)) + qtgf
qsh = csh * (logng(ix,iy,igsp)-logng(ix2,iy,igsp)) + qtgf
qr = abs(qsh/qfl)
c... Because guard-cell values may be distorted from B.C., possibly omit terms on
c... boundary face - shouldnt matter(just set BC) except for guard-cell values
Expand All @@ -6785,7 +6785,7 @@ c physically meaningful gas variables (flux, velocity, etc)
c... now add the convective velocity for charge-exchange neutrals
if(igsp .eq. 1) floxg(ix,iy) =
. floxg(ix,iy) + cngflox(1)*sx(ix,iy)*uu(ix,iy,1)
floxg(ix,iy) = floxg(ix,iy)*2/(lng(ix,iy,igsp)+lng(ix2,iy,igsp))
floxg(ix,iy) = floxg(ix,iy)*2/(logng(ix,iy,igsp)+logng(ix2,iy,igsp))

887 continue
conxg(nx+1,iy) = 0
Expand Down Expand Up @@ -6847,7 +6847,7 @@ c physically meaningful gas variables (flux, velocity, etc)
if(methgy.ne.2) nconv =
. ngy0(ix,iy,igsp)*0.5*(1+sign(1.,qtgf)) +
. ngy1(ix,iy,igsp)*0.5*(1-sign(1.,qtgf))
qsh = csh * (lng(ix,iy,igsp)-lng(ix,iy+1,igsp)) + qtgf
qsh = csh * (logng(ix,iy,igsp)-logng(ix,iy+1,igsp)) + qtgf
qr = abs(qsh/qfl)
if(iy.eq.0 .and. iymnbcl.eq.1) then
qr = gcfacgy*qr
Expand All @@ -6866,7 +6866,7 @@ c physically meaningful gas variables (flux, velocity, etc)
if(igsp .eq. 1) floyg(ix,iy) =
. floyg(ix,iy)+cngfloy(1)*sy(ix,iy)*vy(ix,iy,1)

floyg(ix,iy)=floyg(ix,iy)*2/(lng(ix,iy,igsp)+lng(ix,iy+1,igsp))
floyg(ix,iy)=floyg(ix,iy)*2/(logng(ix,iy,igsp)+logng(ix,iy+1,igsp))
889 continue
890 continue

Expand All @@ -6875,7 +6875,7 @@ c physically meaningful gas variables (flux, velocity, etc)
* --------------------------------------------------------------------

call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
. lng(0:nx+1,0:ny+1,igsp),flngx(0:nx+1,0:ny+1,igsp),
. logng(0:nx+1,0:ny+1,igsp),flngx(0:nx+1,0:ny+1,igsp),
. flngy(0:nx+1,0:ny+1,igsp),0,methg)

c... Addition for nonorthogonal mesh
Expand All @@ -6896,16 +6896,16 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
if ( (ix==ixlb(jx).and.ixmnbcl==1) .or.
. (ix==ixrb(jx).and.ixmxbcl==1) )isxyfl = .false.
enddo
grdnv =( (fym (ix,iy,1)*lng(ix2,iy1 ,igsp) +
. fy0 (ix,iy,1)*lng(ix2,iy ,igsp) +
. fyp (ix,iy,1)*lng(ix2,iy+1,igsp) +
. fymx(ix,iy,1)*lng(ix ,iy1 ,igsp) +
. fypx(ix,iy,1)*lng(ix, iy+1,igsp))
. - (fym (ix,iy,0)*lng(ix ,iy1 ,igsp) +
. fy0 (ix,iy,0)*lng(ix ,iy ,igsp) +
. fyp (ix,iy,0)*lng(ix ,iy+1,igsp) +
. fymx(ix,iy,0)*lng(ix4,iy1 ,igsp) +
. fypx(ix,iy,0)*lng(ix6,iy+1,igsp)) ) /
grdnv =( (fym (ix,iy,1)*logng(ix2,iy1 ,igsp) +
. fy0 (ix,iy,1)*logng(ix2,iy ,igsp) +
. fyp (ix,iy,1)*logng(ix2,iy+1,igsp) +
. fymx(ix,iy,1)*logng(ix ,iy1 ,igsp) +
. fypx(ix,iy,1)*logng(ix, iy+1,igsp))
. - (fym (ix,iy,0)*logng(ix ,iy1 ,igsp) +
. fy0 (ix,iy,0)*logng(ix ,iy ,igsp) +
. fyp (ix,iy,0)*logng(ix ,iy+1,igsp) +
. fymx(ix,iy,0)*logng(ix4,iy1 ,igsp) +
. fypx(ix,iy,0)*logng(ix6,iy+1,igsp)) ) /
. dxnog(ix,iy)

difgx2 = ave( tg(ix ,iy,igsp)/nuix(ix ,iy,igsp),
Expand All @@ -6914,7 +6914,7 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
. 0.5*(nuiz(ix,iy,igsp)+nuiz(ix2,iy,igsp))

flngxy(ix,iy,igsp) = difgx2*(grdnv/cosangfx(ix,iy) -
. (lng(ix2,iy,igsp) - lng(ix,iy,igsp))*
. (logng(ix2,iy,igsp) - logng(ix,iy,igsp))*
. gxf(ix,iy) ) * sx(ix,iy)

c... Now flux limit with flalfgxy
Expand Down Expand Up @@ -6977,14 +6977,14 @@ cc uu(ix,iy,iigsp) = uug(ix,iy,igsp)
do ix = i1, i5
ix2 = ixp1(ix,iy)
fngx(ix,iy,igsp) = flngx(ix,iy,igsp) *
. exp(0.5*(lng(ix,iy,igsp)+lng(ix2,iy,igsp)))
. exp(0.5*(logng(ix,iy,igsp)+logng(ix2,iy,igsp)))
enddo
enddo
c ... now do fniy
do iy = j1, j5 # same loop ranges as for fngy in fd2tra
do ix = i4, i8
fngy(ix,iy,igsp) = flngy(ix,iy,igsp) *
. exp(0.5*(lng(ix,iy,igsp)+lng(ix,iy+1,igsp)))
. exp(0.5*(logng(ix,iy,igsp)+logng(ix,iy+1,igsp)))
enddo
enddo

Expand Down Expand Up @@ -7340,16 +7340,16 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
ix5 = ixm1(ix,iy+1)
ix6 = ixp1(ix,iy+1)
if (methgx .eq. 6) then # log interpolation
grdnv =( exp(fym (ix,iy,1)*log(ng(ix2,iy1 ,igsp)) +
. fy0 (ix,iy,1)*log(ng(ix2,iy ,igsp)) +
. fyp (ix,iy,1)*log(ng(ix2,iy+1,igsp)) +
. fymx(ix,iy,1)*log(ng(ix ,iy1 ,igsp)) +
. fypx(ix,iy,1)*log(ng(ix, iy+1,igsp)))
. - exp(fym (ix,iy,0)*log(ng(ix ,iy1 ,igsp)) +
. fy0 (ix,iy,0)*log(ng(ix ,iy ,igsp)) +
. fyp (ix,iy,0)*log(ng(ix ,iy+1,igsp)) +
. fymx(ix,iy,0)*log(ng(ix4,iy1 ,igsp)) +
. fypx(ix,iy,0)*log(ng(ix6,iy+1,igsp))) ) /
grdnv =( exp(fym (ix,iy,1)*logng(ix2,iy1 ,igsp) +
. fy0 (ix,iy,1)*logng(ix2,iy ,igsp) +
. fyp (ix,iy,1)*logng(ix2,iy+1,igsp) +
. fymx(ix,iy,1)*logng(ix ,iy1 ,igsp) +
. fypx(ix,iy,1)*logng(ix, iy+1,igsp) )
. - exp(fym (ix,iy,0)*logng(ix ,iy1 ,igsp) +
. fy0 (ix,iy,0)*logng(ix ,iy ,igsp) +
. fyp (ix,iy,0)*logng(ix ,iy+1,igsp) +
. fymx(ix,iy,0)*logng(ix4,iy1 ,igsp) +
. fypx(ix,iy,0)*logng(ix6,iy+1,igsp)) ) /
. dxnog(ix,iy)
elseif (methgx .eq. 7) then # inverse interpolation
grdnv =( 1/(fym (ix,iy,1)/ng(ix2,iy1 ,igsp) +
Expand Down
3 changes: 1 addition & 2 deletions bbb/odesetup.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1608,7 +1608,6 @@ call remark('*** Error: ngs <= 0 ***')
call xerrab("")
endif
logng(ix,iy,igsp)=LOG(ABS(ng(ix,iy,igsp)))
lng(ix,iy,igsp) = LOG(ng(ix,iy,igsp))
tg(ix,iy,igsp) = tgs(ix,iy,igsp)
logtg(ix,iy,igsp) = LOG(tg(ix,iy,igsp))
enddo
Expand Down Expand Up @@ -1777,7 +1776,7 @@ call intpvar (tgs(0:,0:,igsp), tg(0:,0:,igsp), 0, nxold, nyold)
ng(ix,iy,igsp) = max(ng(ix,iy,igsp),
. 1.0e-01*ngbackg(igsp))
if (ineudif .eq. 2) then
lng(ix,iy,igsp) = log(ng(ix,iy,igsp))
logng(ix,iy,igsp) = log(ng(ix,iy,igsp))
endif
endif
enddo
Expand Down
4 changes: 2 additions & 2 deletions ppp/debug_parallel.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -443,8 +443,8 @@ subroutine DebugHelper(FileName)
call WriteArrayReal(kye_use,size(kye_use),iunit)
write(iunit,*) "kyi_use"
call WriteArrayReal(kyi_use,size(kyi_use),iunit)
write(iunit,*) "lng"
call WriteArrayReal(lng,size(lng),iunit)
write(iunit,*) "logng"
call WriteArrayReal(logng,size(logng),iunit)
write(iunit,*) "loglambda"
call WriteArrayReal(loglambda,size(loglambda),iunit)
write(iunit,*) "msor"
Expand Down

0 comments on commit ee50085

Please sign in to comment.