Skip to content

Commit

Permalink
Adding evaporation terms for the limiter per T.Rognlien's email of 5-…
Browse files Browse the repository at this point in the history 
…jan-2024
  • Loading branch information
@umansky
umansky committed Jan 6, 2024
1 parent b1b9f99 commit abc7d61
Show file tree
Hide file tree
Showing 2 changed files with 64 additions and 33 deletions.
36 changes: 22 additions & 14 deletions bbb/bbb.v
View file
Original file line number Diff line number Diff line change
Expand Up @@ -331,8 +331,8 @@ cfvgpgy(1:ngspmx) real /ngspmx*0./ +input # Coefs for v*grad(pg) term (poloidal
cfbgt real /0./ +input #Coef for the B x Grad(T) terms.
cfjhf real /1./ +input #Coef for convective cur (fqp) heat flow
jhswitch integer /0/ +input #Coef for the Joule-heating terms
oldseec integer /1/ +input #Switch for Joule-heating bugfix
override integer /0/ +input #Switch to manually override checks on old models
oldseec integer /0/ +input #Switch for Joule-heating bugfix
override integer /1/ +input #Switch to manually override checks on old models
cf2ef real /0./ +input #Coef for ExB drift in 2-direction
cfyef real /0./ +input #Coef for ExB drift in y-direction
cftef real /0./ +input #Coef for ExB drift in toroidal direction
Expand Down Expand Up @@ -899,12 +899,12 @@ albrb(0:ny+1,ngspmx,nxptmx) _real +input #outer plate albedo; used if <1
albedo_by_user integer /0/ +input #if=1, user fills albedoo,i & albdlb,rb
fngxslb(0:ny+1,ngspmx,nxptmx) _real [1/s]+input #inner plt liq vapor gas sour. if sputtlb>0
fngxsrb(0:ny+1,ngspmx,nxptmx) _real [1/s]+input #outer plt liq vapor gas sour. if sputtlb>0
fngxsllim(0:ny+1,ngspmx) _real [1/s] +input #left limiter liq vapor gas sour. if sputtlb>0
fngxsrlim(0:ny+1,ngspmx) _real [1/s] +input #right limiter liq vapor gas sour. if sputtlb>0
fngxsllim(0:ny+1,ngspmx) _real [1/s]+input #left limiter liq vapor gas sour. if sputtlb>0
fngxsrlim(0:ny+1,ngspmx) _real [1/s]+input #right limiter liq vapor gas sour. if sputtlb>0
fngxlb_use(0:ny+1,ngspmx,nxptmx) _real [1/s] +input #user external left plate source
fngxrb_use(0:ny+1,ngspmx,nxptmx) _real [1/s] +input #user external left plate source
fngxllim_use(0:ny+1,ngspmx) _real [1/s] +input #user external left limiter source
fngxrlim_use(0:ny+1,ngspmx) _real [1/s] +input #user external right limiter source
fngxllim_use(0:ny+1,ngspmx) _real [1/s] +input #user external left limiter source
fngxrlim_use(0:ny+1,ngspmx) _real [1/s] +input #user external right limiter source
adatlb(ngspmx,50,nxptmx) _real /1./ +maybeinput #inner albdedo data for each ydati
adatrb(ngspmx,50,nxptmx) _real /1./ +maybeinput #outer albdedo data for each ydati
recycw(ngspmx) real /ngspmx*1e-10/ +input #recycling coef. at side walls
Expand Down Expand Up @@ -944,16 +944,24 @@ ygaslb(10,nxptmx) _real [m] /0./ +input #loc of left-plate sources wrt str
ygasrb(10,nxptmx) _real [m] /0./ +input #loc of right-plate sources wrt strike pt
wgaslb(10,nxptmx) _real [m] /100./ +input #total cos width of left-plate gas sources
wgasrb(10,nxptmx) _real [m] /100./ +input #total cos width of right-plate gas sources
fvaplb(ngspmx,nxptmx) _real /0./ +input #scale factor left-plate evap vapor source
fvaplb(ngspmx,nxptmx) _real /0./ +input #scale factor left-plate evap vapor source
avaplb(ngspmx,nxptmx) _real [k**.5/(m**2s)] /1./ +input #lin coeff left-plate evapor sor
bvaplb(ngspmx,nxptmx) _real [K] /1./ +input #expon. coeff. left-plate evap vapor source
fvaprb(ngspmx,nxptmx) _real /0./ +input #scale factor right-plate evap vapor source
bvaplb(ngspmx,nxptmx) _real [K] /1./ +input #expon. coeff. left-plate evap vapor source
fvaprb(ngspmx,nxptmx) _real /0./ +input #scale factor right-plate evap vapor source
avaprb(ngspmx,nxptmx) _real [k**.5/(m**2s)] /1./ +input #lin coeff right-plate evapor sor
bvaprb(ngspmx,nxptmx) _real [K] /1./ +input #expon coeff. right-plate evap vapor source
tvaplb(0:ny+1,nxptmx) _real [K] +input #left-plate temp for evap; input after alloc
tvaprb(0:ny+1,nxptmx) _real [K] +input #right-plate temp for evap; input after alloc
tvliml(0:ny+1) _real /300./ [K] +input #user left limiter face temp
tvlimr(0:ny+1) _real /300./ [K] +input #user right limiter face temp
bvaprb(ngspmx,nxptmx) _real [K] /1./ +input #expon coeff. right-plate evap vapor source
tvaplb(0:ny+1,nxptmx) _real [K] +input #left-plate temp for evap; input after alloc
tvaprb(0:ny+1,nxptmx) _real [K] +input #right-plate temp for evap; input after alloc
fvapllim(ngspmx) _real /0./ +input #scale factor left-plate evap vapor source
avapllim(ngspmx) _real [k**.5/(m**2s)] /1./ +input #lin coeff left-plate evapor sor
bvapllim(ngspmx) _real [K] /1./ +input #expon. coeff. left-plate evap vapor source
fvaprlim(ngspmx) _real /0./ +input #scale factor right-plate evap vapor source
avaprlim(ngspmx) _real [k**.5/(m**2s)] /1./ +input #lin coeff right-plate evapor sor
bvaprlim(ngspmx) _real [K] /1./ +input #expon coeff. right-plate evap vapor source
tvapllim(0:ny+1) _real [K] +input #left-lim-face temp for evap; input after alloc
tvaprlim(0:ny+1) _real [K] +input #right-lim-fac temp for evap; input after alloc
tvliml(0:ny+1) _real /300./ [K] +input #user left limiter face temp
tvlimr(0:ny+1) _real /300./ [K] +input #user right limiter face temp
isextpltmod integer /0/ +input #=1 use ext gas plate fluxes fngxextlb,rb
# and feixextlb,rb
isextwallmod integer /0/ +input #=1 use ext gas wall fluxes fngyexti,o
Expand Down
61 changes: 42 additions & 19 deletions bbb/boundary.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1655,7 +1655,7 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0
yldot(iv) = -nurlxg * ( fngx(0,iy,igsp) -
. fngxlb_use(iy,igsp,1) +
. fngxslb(iy,igsp,1) + recylb(iy,igsp,1)*flux_inc +
. (1-alblb(iy,igsp,1))*ng(1,iy,igsp)*vxn*sx(0,iy) )
. (1-alblb(iy,igsp,1))*ng(1,iy,igsp)*vxn*sxnp(0,iy) )
. / (vpnorm*n0g(igsp)*sx(0,iy))
endif
if (is1D_gbx.eq.1) yldot(iv) = nurlxg*(ng(1,iy,igsp) -
Expand Down Expand Up @@ -1728,7 +1728,7 @@ ccc yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0
yldot(iv1) = -nurlxg *
. (fnix(ixt,iy,ifld) + recylb(iy,1,jx)*fnix(ixt,iy,1) -
. fngxlb_use(iy,1,jx) +
. (1-alblb(iy,1,jx))*ni(ixt1,iy,ifld)*vxn*sx(ixt,iy) -
. (1-alblb(iy,1,jx))*ni(ixt1,iy,ifld)*vxn*sxnp(ixt,iy) -
. fngxslb(iy,1,jx) ) / (vpnorm*n0(ifld)*sx(ixt,iy))
elseif (recylb(iy,1,jx) <= 0. .and.
. recylb(iy,1,jx) >= -1.) then # recylb is albedo
Expand Down Expand Up @@ -1999,7 +1999,7 @@ cc if (recyce .le. 0) bcen = 0. # gets back to old case
yldot(iv) = -nurlxg * ( fngx(ixt,iy,igsp) -
. fngxlb_use(iy,igsp,jx) -
. fngxslb(iy,igsp,jx) + recylb(iy,igsp,jx)*flux_inc +
. (1-alblb(iy,igsp,jx))*ng(ixt1,iy,igsp)*vxn*sx(ixt,iy) )
. (1-alblb(iy,igsp,jx))*ng(ixt1,iy,igsp)*vxn*sxnp(ixt,iy) )
. / (vpnorm*n0g(igsp)*sx(ixt,iy))
elseif (recylb(iy,igsp,jx) <= 0. .and.
. recylb(iy,igsp,jx) >= -1.) then # recylb is albedo
Expand Down Expand Up @@ -2085,7 +2085,7 @@ call sputchem (isch_sput(igsp),t0p/ev,tvplatlb(iy,jx),
zflux = - sputtlb(iy,igsp,jx) * hflux -
. sputflxlb(iy,igsp,jx) -
. recylb(iy,igsp,jx) * zflux -
. (1-alblb(iy,igsp,jx))*ng(ixt1,iy,igsp)*vxn*sx(ixt1,iy)-
. (1-alblb(iy,igsp,jx))*ng(ixt1,iy,igsp)*vxn*sxnp(ixt1,iy)-
. zflux_chm + fngxslb(iy,igsp,jx)+fngxlb_use(iy,igsp,jx)
yldot(iv) = -nurlxg * (fngx(ixt,iy,igsp) - zflux) /
. (n0(igsp) * vpnorm * sx(ixt,iy))
Expand Down Expand Up @@ -2274,7 +2274,7 @@ ccc yldot(iv) = -nurlxp*(phi(0,iy)-phi(1,iy))/temp0
yldot(iv) = -nurlxg * ( fngx(nx,iy,igsp) +
. fngxrb_use(iy,igsp,1) -
. fngxsrb(iy,igsp,1) + recyrb(iy,igsp,1)*flux_inc -
. (1-albrb(iy,igsp,nxpt))*ng(nx,iy,igsp)*vxn*sx(nx,iy) )
. (1-albrb(iy,igsp,nxpt))*ng(nx,iy,igsp)*vxn*sxnp(nx,iy) )
. / (vpnorm*n0g(igsp)*sx(nx,iy))
endif
endif
Expand Down Expand Up @@ -2347,7 +2347,7 @@ ccc yldot(iv) = -nurlxp*(phi(0,iy)-phi(1,iy))/temp0
yldot(iv1) = nurlxg *
. (fnix(ixt1,iy,ifld) + recyrb(iy,1,jx)*fnix(ixt1,iy,1) +
. fngxrb_use(iy,1,jx) -
. (1-albrb(iy,1,jx))*ni(ixt1,iy,ifld)*vxn*sx(ixt1,iy)
. (1-albrb(iy,1,jx))*ni(ixt1,iy,ifld)*vxn*sxnp(ixt1,iy)
. - fngxsrb(iy,1,jx) ) / (vpnorm*n0(ifld)*sx(ixt1,iy))
elseif (recyrb(iy,1,jx) <= 0. .and.
. recyrb(iy,1,jx) >= -1.) then # recyrb is albedo
Expand Down Expand Up @@ -2636,9 +2636,9 @@ cc if (recyce .le. 0) bcen = 0. # gets back to old case
t0 = max(tg(ixt1,iy,igsp), temin*ev)
vxn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) )
yldot(iv) = nurlxg * ( fngx(ixt1,iy,igsp) +
. fngxrb_use(iy,igsp,jx) -
. fngxrb_use(iy,igsp,jx) -
. fngxsrb(iy,igsp,jx) + recyrb(iy,igsp,jx)*flux_inc -
. (1-albrb(iy,igsp,jx))*ng(ixt1,iy,igsp)*vxn*sx(ixt1,iy) )
. (1-albrb(iy,igsp,jx))*ng(ixt1,iy,igsp)*vxn*sxnp(ixt1,iy) )
. / (vpnorm*n0g(igsp)*sx(ixt1,iy))
elseif (recyrb(iy,igsp,jx) <= 0. .and.
. recyrb(iy,igsp,jx) >= -1.) then # recyrb is albedo
Expand Down Expand Up @@ -2722,7 +2722,7 @@ call sputchem (isch_sput(igsp),t0p/ev,tvplatrb(iy,jx),
zflux = - sputtrb(iy,igsp,jx) * hflux -
. sputflxrb(iy,igsp,jx) -
. recyrb(iy,igsp,jx) * zflux +
. (1-albrb(iy,igsp,jx))*ng(ixt1,iy,igsp)*vxn*sx(ixt1,iy)-
. (1-albrb(iy,igsp,jx))*ng(ixt1,iy,igsp)*vxn*sxnp(ixt1,iy)-
. zflux_chm + fngxsrb(iy,igsp,jx)-fngxrb_use(iy,igsp,jx)
yldot(iv) = nurlxg * (fngx(ixt1,iy,igsp) - zflux) /
. (n0(igsp) * vpnorm * sx(ixt1,iy))
Expand Down Expand Up @@ -3850,7 +3850,7 @@ c This subroutine sets components of the array iseqalg(neq) to be unity
* -- local arrays --
real sycosi(10), sycoso(10)

c Calculate distances along left and right boundaries --
c Calculate distances along left and right poloidal boundaries --
do jx = 1, nxpt
yylb(iysptrx1(jx),jx) = - 0.5/
. ( gy(ixlb(jx),iysptrx1(jx)) * cos(vtag(ixlb(jx),iysptrx1(jx))) )
Expand Down Expand Up @@ -3918,7 +3918,7 @@ real sycosi(10), sycoso(10)
fwsoro(ix,isor) = 0.
286 continue
29 continue

c... We allow for 10 separate sources each on the inner and outer wall
c... If nwsor > 10, arrays must be enlarged
if(nwsor .gt. 10) then
Expand Down Expand Up @@ -4174,7 +4174,7 @@ subroutine wsmodo(ig)
implicit none

Use(Dim) # nx,ny,ngsp,nxpt
## Use(Share) # nyomitmx
Use(Share) # nyomitmx,ix_lim,iy_lims
Use(Xpoint_indices) # ixlb,ixrb
## Use(Math_problem_size) # neqmx(for arrays not used here)
## Use(Selec) # ixp1
Expand Down Expand Up @@ -4240,14 +4240,10 @@ real fsorlb(0:ny+1,npltsor), fsorrb(0:ny+1,npltsor)
enddo
enddo

c... For now (092723), set vapor source arrays on limiters to zero
fngxsllim = 0.
fngxsrlim = 0.

c... Allow for 10 separate sources each on the inner and outer plates
c... If npltsor > 10, arrays must be enlarged
if(npltsor .gt. 10) then
call xerrab ('npltsor > 10, must increase wall source arrays')
call xerrab ('npltsor > 10, increase max wall source arrays')
endif

do jx = 1, nxpt
Expand Down Expand Up @@ -4298,14 +4294,41 @@ call xerrab ('npltsor > 10, must increase wall source arrays')
call remark('**ERR: tvaplb or tvaprb = 0; must set positive')
endif
fngxslb(iy,igsp,jx) = fngxslb(iy,igsp,jx) + fvaplb(igsp,jx)*
. sx(ixlb(jx),iy)*avaplb(igsp,jx)*exp(-bvaplb(igsp,jx)/
. sxnp(ixlb(jx),iy)*avaplb(igsp,jx)*exp(-bvaplb(igsp,jx)/
. tvaplb(iy,jx))/sqrt(tvaplb(iy,jx))
fngxsrb(iy,igsp,jx) = fngxsrb(iy,igsp,jx) - fvaprb(igsp,jx)*
. sx(ixrb(jx),iy)*avaprb(igsp,jx)*exp(-bvaprb(igsp,jx)/
. sxnp(ixrb(jx),iy)*avaprb(igsp,jx)*exp(-bvaprb(igsp,jx)/
. tvaprb(iy,jx))/sqrt(tvaprb(iy,jx))
enddo
endif
enddo
enddo

c... Initial limiter evaporation temp & flux arrays; single-limiter only
do iy = 0, ny+1
if(tvaprlim(iy) == 0) tvaprlim(iy)=1. #prevent rate eval outside range
if(tvapllim(iy) == 0) tvapllim(iy)=1.
do igsp = 1, ngsp
fngxsllim(iy,igsp) = 0.
fngxsrlim(iy,igsp) = 0.
enddo
enddo

c ... Compute and add evaporative gas flux from liquid-wall limiter
c ... Note below fvprlim,llim is not added recursively as no limiter
c ... gas sources unlike fngxslb,rb.

do igsp = 1, ngsp
if (fvapllim(igsp) + fvaprlim(igsp) > 1.e-20) then
do iy = iy_lims, ny
fngxsllim(iy,igsp) = -fvapllim(igsp)*
. sxnp(ix_lim-1,iy)*avapllim(igsp)*exp(-bvapllim(igsp)/
. tvapllim(iy))/sqrt(tvapllim(iy))
fngxsrlim(iy,igsp) = fvaprlim(igsp)*
. sxnp(ix_lim+1,iy)*avaprlim(igsp)*exp(-bvaprlim(igsp)/
. tvaprlim(iy))/sqrt(tvaprlim(iy))
enddo
endif
enddo

return
Expand Down

0 comments on commit abc7d61

Please sign in to comment.