From ee50085afcadc6f83a4bf9129bb30a7b6f2567ce Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Tue, 14 Jan 2025 21:06:16 -0800 Subject: [PATCH] Uses precomputed logng - Replaces lng --- bbb/bbb.v | 1 - bbb/convert.m | 20 +++++++---- bbb/oderhs.m | 80 +++++++++++++++++++++--------------------- bbb/odesetup.m | 3 +- ppp/debug_parallel.F90 | 4 +-- 5 files changed, 57 insertions(+), 51 deletions(-) diff --git a/bbb/bbb.v b/bbb/bbb.v index 0907e6bc..f3047bca 100644 --- a/bbb/bbb.v +++ b/bbb/bbb.v @@ -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; diff --git a/bbb/convert.m b/bbb/convert.m index 1dd8f7ed..b684b97e 100755 --- a/bbb/convert.m +++ b/bbb/convert.m @@ -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)) @@ -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) @@ -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 @@ -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) @@ -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 @@ -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) @@ -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 diff --git a/bbb/oderhs.m b/bbb/oderhs.m index 4c188454..392323fa 100755 --- a/bbb/oderhs.m +++ b/bbb/oderhs.m @@ -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) + @@ -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) - @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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), @@ -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 @@ -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 @@ -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) + diff --git a/bbb/odesetup.m b/bbb/odesetup.m index 8f5bff66..b8947c0a 100755 --- a/bbb/odesetup.m +++ b/bbb/odesetup.m @@ -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 @@ -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 diff --git a/ppp/debug_parallel.F90 b/ppp/debug_parallel.F90 index 82a647d1..4f32443b 100644 --- a/ppp/debug_parallel.F90 +++ b/ppp/debug_parallel.F90 @@ -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"