diff --git a/bbb/bbb.v b/bbb/bbb.v index ca1231a3..eb14f527 100644 --- a/bbb/bbb.v +++ b/bbb/bbb.v @@ -861,62 +861,67 @@ cftgtipltr(ngspmx) real /ngspmx*1./ +threadprivate #right plate Tg B.C.: t cgengmpl real /1./ +input #scale fac mol plate eng loss for Maxw cgengmw real /1./ +input #scale fac mol wall eng loss for Maxw isvacuummodel(ngspmx) integer /6*0/ +input #user:(=0) default model for neutral outer boundary conditions - # (=1) using DEGAS2-calculated tele-transport matrix, half-Mxwellian for outgoing neutrals - # (=2) CX-diffusion for outgoing neutrals - # (=3) combination for outgoing neutrals -cftelematrix(0:nx+1,0:nx+1,ngspmx) _real /0./ +maybeinput # tele-transport matrix if isvacuummodel > 0 -cfpuffmatrix(0:nx+1,ngspmx) _real /0./ +maybeinput # placeholder for puffing matrix if isvacuummodel > 0 - # assuming point puff source, so, the matrix shape is (1,nx) -cfteleout real /1./ +input # scale fac tele-transport neutral source for half-Maxw -fngyteleout(0:nx+1,ngspmx) _real +threadprivate # flux of tele-transport neutral going out of the outer boundary, -fngytelein(0:nx+1,ngspmx) _real +threadprivate # flux of tele-transport neutral coming back into the outer boundary, + # (=1) Vacuum Neutral Model using tele-transport matrix, assuming Max-wellian out-going neutrals + # (=2) Vacuum Neutral Model using tele-transport matrix, assuming CX-diffusion out-going neutrals + # (=3) Vacuum Neutral Model using tele-transport matrix, combining Max-wellian and CX-diffusion + # for out-going neutrals based on local Knudsen Number. +cftelematrix(1:2,0:nx+1,0:nx+1,ngspmx) _real /0./ +maybeinput # tele-transport matrix if isvacuummodel = 1,2,3 + # (1,:,:,:) for inner PFR boundary + # (2,:,:,:) for outer wall boundary +cfpuffmatrix(1:2,0:nx+1,ngspmx) _real /0./ +maybeinput # placeholder for puffing matrix if isvacuummodel > 0 + # assuming point puff source, so, the matrix shape is (:,1:nx,:) +cfteleout(1:2) real /2*1./ +input # scale fac tele-transport neutral source for half-Maxw +fngyteleout(1:2,0:nx+1,ngspmx) _real +threadprivate # flux of tele-transport neutral going out of the outer boundary, +fngytelein(1:2,0:nx+1,ngspmx) _real +threadprivate # flux of tele-transport neutral coming back into the outer boundary, # = fngyteleout*cftelematrix -fngytelemaxw(0:nx+1,ngspmx) _real +threadprivate # flux of tele-transport neutral going out of the outer boundary, assuming half-Maxw -fngytelediff(0:nx+1,ngspmx) _real +threadprivate # flux of tele-transport neutral going out of the outer boundary, assuming CX diffusion +fngytelemaxw(1:2,0:nx+1,ngspmx) _real +threadprivate # flux of tele-transport neutral going out of the outer boundary, assuming half-Maxw + # positive sign means that flux is pointing out of the boundary for both i,o boundaries. +fngytelediff(1:2,0:nx+1,ngspmx) _real +threadprivate # flux of tele-transport neutral going out of the outer boundary, assuming CX diffusion # Three compoments: fngytelediff1 + fngytelediff2 + fngytelediff3 # 1) fngytelediff1 = c1*ng # 2) fngytelediff2 = c2*dng/dtau, where dng/dtau is gradient of ng parallel to the wall # 3) fngytelediff3 = c3*dng/dnu, where dng/dnu is gradient of ng along the wall normal -fngytelediff1(0:nx+1,ngspmx) _real +threadprivate # first component of fngytelediff -fngytelediff2(0:nx+1,ngspmx) _real +threadprivate # second component of fngytelediff -fngytelediff3(0:nx+1,ngspmx) _real +threadprivate # third component of fngytelediff -fngytelediff30(0:nx+1,ngspmx) _real +threadprivate # third component of fngytelediff before applying limiter + # positive sign means that flux is pointing out of the boundary for both i,o boundaries. +fngytelediff1(1:2,0:nx+1,ngspmx) _real +threadprivate # first component of fngytelediff +fngytelediff2(1:2,0:nx+1,ngspmx) _real +threadprivate # second component of fngytelediff +fngytelediff3(1:2,0:nx+1,ngspmx) _real +threadprivate # third component of fngytelediff +fngytelediff30(1:2,0:nx+1,ngspmx) _real +threadprivate # third component of fngytelediff before applying limiter cftelelim real [m] /1./ +input # scale factor for fngytelediff3 limiter, limiting the maximum fngytelediff3 below lambda_i/lambda_n ~ 1.0 ntelelim real /2./ +input # exponent for fngytelediff3 limiter lteleout real [m] /1./ +input # scale length used to calculate neutral Knusen number, Kn = lmfp/lteleout -Kn(0:nx+1,ngspmx) _real +threadprivate # Knusen number at the outer boundary, = lmfpn/lteleout, where lmfpn = sqrt(Tg/mg)/nucx +Kn(1:2,0:nx+1,ngspmx) _real +threadprivate # Knusen number at the outer boundary, = lmfpn/lteleout, where lmfpn = sqrt(Tg/mg)/nucx Knb1 real /0.01/ +input # lower bound of Knusen number Knb2 real /1./ +input # upper bound of Knusen number # if Kn <= Knb1, fngyteleout = fngytelediff # else if Kn >= Knb2, fngyteleout = fngytelemaxw # else fngyteleout = c1*fngytelemaxw + c2*fngytelediff # where c1 = (Kn-Knb1)/(Knb2-Knb1), and c2 = 1-c1 -fmgyteleout(0:nx+1,ngsp) _real +threadprivate # mom. flux of tele-transport neutral going out of the outer boundary, -fmgytelein(0:nx+1,ngsp) _real +threadprivate # mom. flux of tele-transport neutral coming back into the outer boundary, +fmgyteleout(1:2,0:nx+1,ngsp) _real +threadprivate # mom. flux of tele-transport neutral going out of the outer boundary, +fmgytelein(1:2,0:nx+1,ngsp) _real +threadprivate # mom. flux of tele-transport neutral coming back into the outer boundary, # = 0. for current model, assuming wall absorb all the mom. and energy -fmgytelemaxw(0:nx+1,ngsp) _real +threadprivate # mom. flux of tele-transport neutral going out of the outer boundary, assuming half-Maxw +fmgytelemaxw(1:2,0:nx+1,ngsp) _real +threadprivate # mom. flux of tele-transport neutral going out of the outer boundary, assuming half-Maxw # fegytelemaxw = fngytelemaxw*mu_p -fmgytelediff(0:nx+1,ngsp) _real +threadprivate # mom. flux of tele-transport neutral going out of the outer boundary, assuming CX diffusion +fmgytelediff(1:2,0:nx+1,ngsp) _real +threadprivate # mom. flux of tele-transport neutral going out of the outer boundary, assuming CX diffusion # Three compoments: fmgytelediff1 + fmgytelediff2 + fmgytelediff3 # 1) fmgytelediff1 = fngytelediff1*mu_p # 2) fmgytelediff2 = fngytelediff2*mu_p + 'something' ('something' ignored for current model) # 3) fmgytelediff3 = fngytelediff3*mu_p -fmgytelediff1(0:nx+1,ngsp) _real +threadprivate # first component of fmgytelediff -fmgytelediff2(0:nx+1,ngsp) _real +threadprivate # second component of fmgytelediff -fmgytelediff3(0:nx+1,ngsp) _real +threadprivate # third component of fmgytelediff -fegyteleout(0:nx+1,ngsp) _real +threadprivate # eng. flux of tele-transport neutral going out of the outer boundary, -fegytelein(0:nx+1,ngsp) _real +threadprivate # eng. flux of tele-transport neutral coming back into the outer boundary, +fmgytelediff1(1:2,0:nx+1,ngsp) _real +threadprivate # first component of fmgytelediff +fmgytelediff2(1:2,0:nx+1,ngsp) _real +threadprivate # second component of fmgytelediff +fmgytelediff3(1:2,0:nx+1,ngsp) _real +threadprivate # third component of fmgytelediff +fegyteleout(1:2,0:nx+1,ngsp) _real +threadprivate # eng. flux of tele-transport neutral going out of the outer boundary, +fegytelein(1:2,0:nx+1,ngsp) _real +threadprivate # eng. flux of tele-transport neutral coming back into the outer boundary, # = 0. for current model, assuming wall absorb all the mom. and energy -fegytelemaxw(0:nx+1,ngsp) _real +threadprivate # eng. flux of tele-transport neutral going out of the outer boundary, assuming half-Maxw +fegytelemaxw(1:2,0:nx+1,ngsp) _real +threadprivate # eng. flux of tele-transport neutral going out of the outer boundary, assuming half-Maxw # fegytelemaxw = fngytelemaxw*(0.5mu_p^2 + 2T), omit 0.5mu_p^2 for internal -fegytelediff(0:nx+1,ngsp) _real +threadprivate # eng. flux of tele-transport neutral going out of the outer boundary, assuming CX diffusion +fegytelediff(1:2,0:nx+1,ngsp) _real +threadprivate # eng. flux of tele-transport neutral going out of the outer boundary, assuming CX diffusion # Three compoments: fegytelediff1 + fegytelediff2 + fegytelediff3 # 1) fegytelediff1 = fngytelediff1*(0.5mu_p^2 + 2Ti), omit 0.5mu_p^2 for internal # 2) fegytelediff2 = fngytelediff2*(0.5mu_p^2 + 5/2Ti), omit 0.5mu_p^2 for internal # 3) fegytelediff3 = fngytelediff3*(0.5mu_p^2 + 5/2Ti), omit 0.5mu_p^2 for internal -fegytelediff1(0:nx+1,ngsp) _real +threadprivate # first component of fegytelediff -fegytelediff2(0:nx+1,ngsp) _real +threadprivate # second component of fegytelediff -fegytelediff3(0:nx+1,ngsp) _real +threadprivate # third component of fegytelediff +fegytelediff1(1:2,0:nx+1,ngsp) _real +threadprivate # first component of fegytelediff +fegytelediff2(1:2,0:nx+1,ngsp) _real +threadprivate # second component of fegytelediff +fegytelediff3(1:2,0:nx+1,ngsp) _real +threadprivate # third component of fegytelediff ***** Outpwall: # Arrays used to communicate wall fluxes to a wall simulation code diff --git a/bbb/pandf_boundary.m b/bbb/pandf_boundary.m index 147e82c0..82967aae 100755 --- a/bbb/pandf_boundary.m +++ b/bbb/pandf_boundary.m @@ -1338,6 +1338,7 @@ real yldot(neq) . feiytotc, fmiytotc, fngytotc, fniytotc, fqytotc, sytotc, . sfeeytotc, sfeiytotc, lengg, uztotc, uztotc1, uztotc2, . upbdry1, upbdry2, upbdry, r_major, xtotc, zflux, ut0 + real Mi_nu,Mi_taunu,Mi_nunu,dng2dnu,dng2dtau,lmfpn,fteledifftrim,fngytelelim real ave, t0,t1 ave(t0,t1) = 2*t0*t1 / (cutlo+t0+t1) @@ -1367,6 +1368,66 @@ real yldot(neq) # this is a very large if-loop; check for 'j3 index' to find end do ifld = 1, nisp + if (isupgon(1) .eq. 1 .and. zi(ifld) .eq. 0.0 .and. isvacuummodel(1) .gt. 0) then #..zml + do ix = i4+1-ixmnbcl, i8-1+ixmxbcl + if (isixcore(ix) .ne. 1) then + # half-Maxwellian + t0 = max(tg(ix,1,1),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mi(ifld)) ) + fngytelemaxw(1,ix,1) = cfteleout(1)*ni(ix,1,ifld)*vyn*sy(ix,0) + # CX diffusion + ix1 = ixm1(ix,1) + ix2 = ixp1(ix,1) + t1 = max(ti(ix,1),temin*ev) + Mi_nu = -sqrt( 0.5*t1/(pi*mi(1)) ) + Mi_taunu = 0.5*( up(ix,1,0)*rrv(ix,1) + up(ix1,1,0)*rrv(ix1,1) )*Mi_nu # = Utau*Mi_nu + Mi_nunu = 0.5*t0/mi(1) + dng2dnu = 0.5*( (ngy1(ix,0,1) - ngy0(ix,0,1))/dynog(ix,0) + . + (ngy1(ix,1,1) - ngy0(ix,1,1))/dynog(ix,1) ) + dng2dtau = 0.5*( (ng(ix2,1,1) - ng(ix,1,1))*gxf(ix,1) + . + (ng(ix,1,1) - ng(ix1,1,1))*gxf(ix1,1) ) + fngytelediff1(1,ix,1) = - nucx(ix,1,1)/(nucx(ix,1,1)+nuiz(ix,1,1)) * + . ng(ix,1,1)*Mi_nu * sy(ix,0) + fngytelediff2(1,ix,1) = nucx(ix,1,1)/(nucx(ix,1,1)+nuiz(ix,1,1))**2 * + . dng2dtau*Mi_taunu * sy(ix,0) + fngytelediff30(1,ix,1) = nucx(ix,1,1)/(nucx(ix,1,1)+nuiz(ix,1,1))**2 * + . dng2dnu*Mi_nunu * sy(ix,0) + # The third component is dominant in divertor regions and can cause numerical instability + # We add a flux limiter to this term fngytelediff3, + # where the constant 1.2533 is a result of pi/sqrt(2*pi) + fngytelelim = cftelelim * 1.2533*abs(nucx(ix,1,1)/(nucx(ix,1,1)+nuiz(ix,1,1))*ng(ix,1,1)*Mi_nu * sy(ix,0)) + fngytelediff3(1,ix,1) = fngytelediff30(1,ix,1)/ + . (1 + (abs(fngytelediff30(1,ix,1)/fngytelelim))**ntelelim )**(1./ntelelim) + fngytelediff(1,ix,1) = cfteleout(1) * (fngytelediff1(1,ix,1) + + . fngytelediff2(1,ix,1) + fngytelediff3(1,ix,1)) + fteledifftrim = max(fngytelediff(1,ix,1),0.) # Trim to avoid negative influx + if (isvacuummodel(1) .eq. 1) then # half-Maxwellian + fngyteleout(1,ix,1) = fngytelemaxw(1,ix,1) + elseif (isvacuummodel(1) .eq. 2) then # CX diff, Eq.(3.88) in AFN manual from KU Lueven + fngyteleout(1,ix,1) = fteledifftrim + elseif (isvacuummodel(1) .eq. 3) then + lmfpn = sqrt(t0/mi(ifld))/nucx(ix,1,1) + Kn(1,ix,1) = lmfpn/lteleout + if (Kn(1,ix,1) <= Knb1) then + fngyteleout(1,ix,1) = fteledifftrim + elseif (Kn(1,ix,1) >= Knb2) then + fngyteleout(1,ix,1) = fngytelemaxw(1,ix,1) + else + fngyteleout(1,ix,1) = (Kn(1,ix,1)-Knb1)/(Knb2-Knb1)*fngytelemaxw(1,ix,1) + + . (Knb2-Kn(1,ix,1))/(Knb2-Knb1)*fteledifftrim + endif + else + write (*,*) 'isvacuummodel > 3 not available for atom at PFR, placeholder for future extention' + endif + endif + enddo + do ix = 0, nx+1 + if (isixcore(ix) .ne. 1) then + fngytelein(1,ix,1) = sum(fngyteleout(1,:,1)*cftelematrix(1,:,ix,1)) + endif + enddo + endif + do ix = i4+1-ixmnbcl, i8-1+ixmxbcl # ix-loop over ion dens if (isnionxy(ix,0,ifld)==1) then iv1 = idxn(ix,0,ifld) @@ -1421,6 +1482,7 @@ call xerrab("*** isngcore=2 option unvailable ***") nharmave = 2.*(ni(ix,0,ifld)*ni(ix,1,ifld)) / . (ni(ix,0,ifld)+ni(ix,1,ifld)) fng_alb = (1-albedoi(ix,1))*nharmave*vyn*sy(ix,0) + if (isvacuummodel(1) .gt. 0) fng_alb = fngyteleout(1,ix,1) - fngytelein(1,ix,1) yldot(iv1) = -nurlxg*( fniy(ix,0,ifld) + fng_alb - . fng_chem ) / . (vyn*sy(ix,0)* n0(ifld)) @@ -1621,6 +1683,44 @@ call xerrab ('** isnicore value not valid option **') call xerrab ("*** Illegal setting of isupcore ***") endif + elseif (isupgon(1) .eq. 1 .and. zi(ifld) .eq. 0.0 .and. isvacuummodel(1) .gt. 0) then # VNM + ix1 = ixm1(ix,ny) + # Maxwellian + fmgytelemaxw(1,ix,1) = 0.5*fngytelemaxw(1,ix,1)*mg(1)*(up(ix,1,ifld)+up(ix1,1,ifld)) + # CX diffusion + fmgytelediff1(1,ix,1) = 0.5*fngytelediff1(1,ix,1)*mg(1)*(up(ix,1,ifld)+up(ix1,1,ifld)) + fmgytelediff2(1,ix,1) = 0.5*fngytelediff2(1,ix,1)*mg(1)*(up(ix,1,ifld)+up(ix1,1,ifld)) #..approximation, omitting small terms + fmgytelediff3(1,ix,1) = 0.5*fngytelediff3(1,ix,1)*mg(1)*(up(ix,1,ifld)+up(ix1,1,ifld)) + fmgytelediff(1,ix,1) = cfteleout(1)*(fmgytelediff1(1,ix,1) + + . fmgytelediff2(1,ix,1) + fmgytelediff3(1,ix,1)) + # Trim diffusion to avoid negative influx + if (fngytelediff(1,ix,1) .gt. 0.) then + fteledifftrim = fmgytelediff(1,ix,1) + else + fteledifftrim = 0. + endif + # compute fmgyteleout + if (isvacuummodel(1) .eq. 1) then + fmgyteleout(1,ix,1) = fmgytelemaxw(1,ix,1) + elseif (isvacuummodel(1) .eq. 2) then + fmgyteleout(1,ix,1) = fteledifftrim + elseif (isvacuummodel(1) .eq. 3) then + if (Kn(1,ix,1) <= Knb1) then + fmgyteleout(1,ix,1) = fteledifftrim + elseif (Kn(1,ix,1) >= Knb2) then + fmgyteleout(1,ix,1) = fmgytelemaxw(1,ix,1) + else + fmgyteleout(1,ix,1) = (Kn(1,ix,1)-Knb1)/(Knb2-Knb1)*fmgytelemaxw(1,ix,1) + + . (Knb2-Kn(1,ix,1))/(Knb2-Knb1)*fteledifftrim + endif + else + write (*,*) 'isvacuummodel > 3 not available for upg at PFR, placeholder for future extention' + endif + # place holder for momentu influx fmgytelein + fmgytelein(1,ix,1) = 0. + yldot(iv2) = nurlxu * (fmiy(ix,0,ifld) + (fmgyteleout(1,ix,1)-fmgytelein(1,ix,1)))/ + . (vpnorm*sy(ix,0)*fnorm(ifld)) + elseif (isupwiix(ix,ifld)==1) then # isixcore if-test; PF bdry yldot(iv2) = -nurlxu * fmiy(ix,0,ifld) / . (vpnorm*sy(ix,0)*fnorm(ifld)) @@ -2058,54 +2158,90 @@ call sputchem (isch_sput(igsp), eng_sput, tvwalli(ix), . (temp0*ev) endif else # PF wall - if (istgpfcix(ix,igsp) == 0) then # fixed Tg - yldot(iv) = nurlxg*(tgwall(igsp)*ev-tg(ix,0,igsp))/ + if (isvacuummodel(igsp) .gt. 0) then + fegytelemaxw(1,ix,igsp) = fngytelemaxw(1,ix,igsp)*2.*tg(ix,1,igsp) + fegytelediff1(1,ix,igsp) = fngytelediff1(1,ix,igsp)*2.*tg(ix,1,igsp) + fegytelediff2(1,ix,igsp) = fngytelediff2(1,ix,igsp)*2.5*tg(ix,1,igsp) + fegytelediff3(1,ix,igsp) = fngytelediff3(1,ix,igsp)*2.5*tg(ix,1,igsp) + fegytelediff(1,ix,igsp) = cfteleout(1)*(fegytelediff1(1,ix,igsp) + + . fegytelediff2(1,ix,igsp) + fegytelediff3(1,ix,igsp)) + # Trim diffusion to avoid negative influx + if (fngytelediff(1,ix,igsp) .gt. 0.) then + fteledifftrim = fegytelediff(1,ix,igsp) + else + fteledifftrim = 0. + endif + # compute fegyteleout + if (isvacuummodel(igsp) .eq. 1) then + fegyteleout(1,ix,igsp) = fegytelemaxw(1,ix,igsp) + elseif (isvacuummodel(igsp) .eq. 2) then + fegyteleout(1,ix,igsp) = fteledifftrim + elseif (isvacuummodel(igsp) .eq. 3) then + if (Kn(1,ix,igsp) <= Knb1) then + fegyteleout(1,ix,igsp) = fteledifftrim + elseif (Kn(1,ix,igsp) >= Knb2) then + fegyteleout(1,ix,igsp) = fegytelemaxw(1,ix,igsp) + else + fegyteleout(1,ix,igsp) = (Kn(1,ix,igsp)-Knb1)/(Knb2-Knb1)*fegytelemaxw(1,ix,igsp) + + . (Knb2-Kn(1,ix,igsp))/(Knb2-Knb1)*fteledifftrim + endif + else + write (*,*) 'isvacuummodel > 3 not available for tg at PFR, placeholder for future extention' + endif + # place holder for momentu influx fmgytelein + fegytelein(1,ix,igsp) = 0. + yldot(iv) = nurlxg*( fegy(ix,0,igsp) + (fegyteleout(1,ix,igsp)-fegytelein(1,ix,igsp)) )/ + . (sy(ix,ny)*vpnorm*ennorm) + else + if (istgpfcix(ix,igsp) == 0) then # fixed Tg + yldot(iv) = nurlxg*(tgwall(igsp)*ev-tg(ix,0,igsp))/ . (temp0*ev) - elseif (istgpfcix(ix,igsp) == 1) # extrapolation - tbound = tg(ix,1,igsp) - gyf(ix,1)* - . (tg(ix,2,igsp)-tg(ix,1,igsp))/gyf(ix,0) - tbound = max(tbound, 0.25*tbmin*ev) #tbmin=.1 eV - yldot(iv) = nurlxi *(tbound - tg(ix,0,igsp))/(temp0*ev) - elseif (istgpfcix(ix,igsp) == 2) # specified gradient - yldot(iv) = nurlxi*( (tg(ix,1,igsp) - tg(ix,0,igsp)) - + elseif (istgpfcix(ix,igsp) == 1) # extrapolation + tbound = tg(ix,1,igsp) - gyf(ix,1)* + . (tg(ix,2,igsp)-tg(ix,1,igsp))/gyf(ix,0) + tbound = max(tbound, 0.25*tbmin*ev) #tbmin=.1 eV + yldot(iv) = nurlxi *(tbound - tg(ix,0,igsp))/(temp0*ev) + elseif (istgpfcix(ix,igsp) == 2) # specified gradient + yldot(iv) = nurlxi*( (tg(ix,1,igsp) - tg(ix,0,igsp)) - . 0.5*(tg(ix,1,igsp) + tg(ix,0,igsp))/ . (gyf(ix,0)*lytg(1,igsp)) )/(temp0*ev) - elseif (istgpfcix(ix,igsp) == 3) #Maxwell thermal flux to wall - t0 = max(cdifg(igsp)*tg(ix,1,igsp), temin*ev) - vyn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) - yldot(iv) = -nurlxg*( fegy(ix,0,igsp) + 2*cgengmw* + elseif (istgpfcix(ix,igsp) == 3) #Maxwell thermal flux to wall + t0 = max(cdifg(igsp)*tg(ix,1,igsp), temin*ev) + vyn = 0.25 * sqrt( 8*t0/(pi*mg(igsp)) ) + yldot(iv) = -nurlxg*( fegy(ix,0,igsp) + 2*cgengmw* . ng(ix,1,igsp)*vyn*t0*sy(ix,0) )/ . (sy(ix,0)*vpnorm*ennorm) - elseif (istgpfcix(ix,igsp) == 4) - t0 = max(tg(ix,0,igsp),tgmin*ev) - vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) - nharmave = 2.*(ng(ix,0,igsp)*ng(ix,1,igsp)) / + elseif (istgpfcix(ix,igsp) == 4) + t0 = max(tg(ix,0,igsp),tgmin*ev) + vyn = sqrt( 0.5*t0/(pi*mg(igsp)) ) + nharmave = 2.*(ng(ix,0,igsp)*ng(ix,1,igsp)) / . (ng(ix,0,igsp)+ng(ix,1,igsp)) - fng_alb = (1-albedoi(ix,1))* + fng_alb = (1-albedoi(ix,1))* . nharmave*vyn*sy(ix,0) - fng_chem = 0. #..double check - yldot(iv)=-nurlxg*(fegy(ix,0,igsp) + cfalbedo*fng_alb*t0 + fng_chem = 0. #..double check + yldot(iv)=-nurlxg*(fegy(ix,0,igsp) + cfalbedo*fng_alb*t0 . - 2.*fng_chem*t0) . /(vpnorm*ennorm*sy(ix,0)) - fniy_recy = 0. - if (matwalli(ix) .gt. 0) then - if (recycwit(ix,igsp,1) .gt. 0) then - fniy_recy = recycwit(ix,igsp,1)*fac2sp*fniy(ix,0,1) - if (isrefluxclip==1) fniy_recy=min(fniy_recy,0.) - yldot(iv)=-nurlxg*(fegy(ix,0,igsp) + cfalbedo*fng_alb*t0 + fniy_recy = 0. + if (matwalli(ix) .gt. 0) then + if (recycwit(ix,igsp,1) .gt. 0) then + fniy_recy = recycwit(ix,igsp,1)*fac2sp*fniy(ix,0,1) + if (isrefluxclip==1) fniy_recy=min(fniy_recy,0.) + yldot(iv)=-nurlxg*(fegy(ix,0,igsp) + cfalbedo*fng_alb*t0 . - 2.*fng_chem*t0 . + fniy_recy*(1.-cfdiss) . *cfalbedo*recycwe . *ti(ix,0) ) . /(vpnorm*ennorm*sy(ix,0)) + endif endif - endif - elseif (istgpfc(igsp) == 5) then # tg=ti*cftgtipfc - yldot(iv) = nurlxg*(ti(ix,0)*cftgtipfc(igsp) + elseif (istgpfc(igsp) == 5) then # tg=ti*cftgtipfc + yldot(iv) = nurlxg*(ti(ix,0)*cftgtipfc(igsp) . -tg(ix,0,igsp))/ . (temp0*ev) - elseif (istgpfc(igsp) > 5) - call xerrab("***Input error: invalid istgpfc ***") + elseif (istgpfc(igsp) > 5) + call xerrab("***Input error: invalid istgpfc ***") + endif endif endif @@ -2442,7 +2578,7 @@ real yldot(neq) # half-Maxwellian t0 = max(tg(ix,ny,1),tgmin*ev) vyn = sqrt( 0.5*t0/(pi*mi(ifld)) ) - fngytelemaxw(ix,1) = cfteleout*ni(ix,ny,ifld)*vyn*sy(ix,ny) + fngytelemaxw(2,ix,1) = cfteleout(2)*ni(ix,ny,ifld)*vyn*sy(ix,ny) # CX diffusion ix1 = ixm1(ix,ny) ix2 = ixp1(ix,ny) @@ -2454,44 +2590,44 @@ real yldot(neq) . + (ngy1(ix,ny-1,1) - ngy0(ix,ny-1,1))/dynog(ix,ny-1) ) dng2dtau = -0.5*( (ng(ix2,ny,1) - ng(ix,ny,1))*gxf(ix,ny) . + (ng(ix,ny,1) - ng(ix1,ny,1))*gxf(ix1,ny) ) - fngytelediff1(ix,1) = - nucx(ix,ny,1)/(nucx(ix,ny,1)+nuiz(ix,ny,1)) * + fngytelediff1(2,ix,1) = - nucx(ix,ny,1)/(nucx(ix,ny,1)+nuiz(ix,ny,1)) * . ng(ix,ny,1)*Mi_nu * sy(ix,ny) - fngytelediff2(ix,1) = nucx(ix,ny,1)/(nucx(ix,ny,1)+nuiz(ix,ny,1))**2 * + fngytelediff2(2,ix,1) = nucx(ix,ny,1)/(nucx(ix,ny,1)+nuiz(ix,ny,1))**2 * . dng2dtau*Mi_taunu * sy(ix,ny) - fngytelediff30(ix,1) = nucx(ix,ny,1)/(nucx(ix,ny,1)+nuiz(ix,ny,1))**2 * + fngytelediff30(2,ix,1) = nucx(ix,ny,1)/(nucx(ix,ny,1)+nuiz(ix,ny,1))**2 * . dng2dnu*Mi_nunu * sy(ix,ny) # The third component is dominant in divertor regions and can cause numerical instability # We add a flux limiter to this term fngytelediff3, # where the constant 1.2533 is a result of pi/sqrt(2*pi) fngytelelim = cftelelim * 1.2533*abs(nucx(ix,ny,1)/(nucx(ix,ny,1)+nuiz(ix,ny,1))*ng(ix,ny,1)*Mi_nu * sy(ix,ny)) - fngytelediff3(ix,1) = fngytelediff30(ix,1)/ - . (1 + (abs(fngytelediff30(ix,1)/fngytelelim))**ntelelim )**(1./ntelelim) - fngytelediff(ix,1) = cfteleout * (fngytelediff1(ix,1) + - . fngytelediff2(ix,1) + fngytelediff3(ix,1)) - fteledifftrim = max(fngytelediff(ix,1),0.) # Trim to avoid negative influx + fngytelediff3(2,ix,1) = fngytelediff30(2,ix,1)/ + . (1 + (abs(fngytelediff30(2,ix,1)/fngytelelim))**ntelelim )**(1./ntelelim) + fngytelediff(2,ix,1) = cfteleout(2) * (fngytelediff1(2,ix,1) + + . fngytelediff2(2,ix,1) + fngytelediff3(2,ix,1)) + fteledifftrim = max(fngytelediff(2,ix,1),0.) # Trim to avoid negative influx if (isvacuummodel(1) .eq. 1) then # half-Maxwellian - fngyteleout(ix,1) = fngytelemaxw(ix,1) + fngyteleout(2,ix,1) = fngytelemaxw(2,ix,1) elseif (isvacuummodel(1) .eq. 2) then # CX diff, Eq.(3.88) in AFN manual from KU Lueven - fngyteleout(ix,1) = fteledifftrim + fngyteleout(2,ix,1) = fteledifftrim elseif (isvacuummodel(1) .eq. 3) then lmfpn = sqrt(t0/mi(ifld))/nucx(ix,ny,1) - Kn(ix,1) = lmfpn/lteleout - if (Kn(ix,1) <= Knb1) then - fngyteleout(ix,1) = fteledifftrim - elseif (Kn(ix,1) >= Knb2) then - fngyteleout(ix,1) = fngytelemaxw(ix,1) + Kn(2,ix,1) = lmfpn/lteleout + if (Kn(2,ix,1) <= Knb1) then + fngyteleout(2,ix,1) = fteledifftrim + elseif (Kn(2,ix,1) >= Knb2) then + fngyteleout(2,ix,1) = fngytelemaxw(2,ix,1) else - fngyteleout(ix,1) = (Kn(ix,1)-Knb1)/(Knb2-Knb1)*fngytelemaxw(ix,1) + - . (Knb2-Kn(ix,1))/(Knb2-Knb1)*fteledifftrim + fngyteleout(2,ix,1) = (Kn(2,ix,1)-Knb1)/(Knb2-Knb1)*fngytelemaxw(2,ix,1) + + . (Knb2-Kn(2,ix,1))/(Knb2-Knb1)*fteledifftrim endif else - write (*,*) 'isvacuummodel > 3 not available, placeholder for future extention' + write (*,*) 'isvacuummodel > 3 not available for atom at outer wall, placeholder for future extention' endif enddo do ix = 0, nx+1 #..compute fngytelein t0 = max(tg(ix,ny,1),tgmin*ev) vyn = sqrt( 0.5*t0/(pi*mi(ifld)) ) - fngytelein(ix,1) = sum(fngyteleout(:,1)*cftelematrix(:,ix,1)) + fngytelein(2,ix,1) = sum(fngyteleout(2,:,1)*cftelematrix(2,:,ix,1)) enddo endif do ix = i4+1-ixmnbcl, i8-1+ixmxbcl @@ -2514,7 +2650,7 @@ real yldot(neq) nharmave = 2.*(ni(ix,ny,ifld)*ni(ix,ny+1,ifld)) / , (ni(ix,ny,ifld)+ni(ix,ny+1,ifld)) fng_alb = (1-albedoo(ix,1))*nharmave*vyn*sy(ix,ny) - if (isvacuummodel(1) .gt. 0) fng_alb = fngyteleout(ix,1) - fngytelein(ix,1) + if (isvacuummodel(1) .gt. 0) fng_alb = fngyteleout(2,ix,1) - fngytelein(2,ix,1) yldot(iv1) = nurlxg*( fniy(ix,ny,ifld) - fng_alb + . fng_chem ) / (vyn*sy(ix,ny)* n0(ifld)) @@ -2609,38 +2745,39 @@ real yldot(neq) if (isupgon(1) .eq. 1 .and. zi(ifld) .eq. 0.0 .and. isvacuummodel(1) .gt. 0) then ix1 = ixm1(ix,ny) # Maxwellian - fmgytelemaxw(ix,1) = 0.5*fngytelemaxw(ix,1)*mg(1)*(up(ix,ny,ifld)+up(ix1,ny,ifld)) + fmgytelemaxw(2,ix,1) = 0.5*fngytelemaxw(2,ix,1)*mg(1)*(up(ix,ny,ifld)+up(ix1,ny,ifld)) # CX diffusion - fmgytelediff1(ix,1) = 0.5*fngytelediff1(ix,1)*mg(1)*(up(ix,ny,ifld)+up(ix1,ny,ifld)) - fmgytelediff2(ix,1) = 0.5*fngytelediff2(ix,1)*mg(1)*(up(ix,ny,ifld)+up(ix1,ny,ifld)) #..approximation, omitting small terms - fmgytelediff3(ix,1) = 0.5*fngytelediff3(ix,1)*mg(1)*(up(ix,ny,ifld)+up(ix1,ny,ifld)) - fmgytelediff(ix,1) = cfteleout*(fmgytelediff1(ix,1) + fmgytelediff2(ix,1) + fmgytelediff3(ix,1)) + fmgytelediff1(2,ix,1) = 0.5*fngytelediff1(2,ix,1)*mg(1)*(up(ix,ny,ifld)+up(ix1,ny,ifld)) + fmgytelediff2(2,ix,1) = 0.5*fngytelediff2(2,ix,1)*mg(1)*(up(ix,ny,ifld)+up(ix1,ny,ifld)) #..approximation, omitting small terms + fmgytelediff3(2,ix,1) = 0.5*fngytelediff3(2,ix,1)*mg(1)*(up(ix,ny,ifld)+up(ix1,ny,ifld)) + fmgytelediff(2,ix,1) = cfteleout(2)*(fmgytelediff1(2,ix,1) + + . fmgytelediff2(2,ix,1) + fmgytelediff3(2,ix,1)) # Trim diffusion to avoid negative influx - if (fngytelediff(ix,1) .gt. 0.) then - fteledifftrim = fmgytelediff(ix,1) + if (fngytelediff(2,ix,1) .gt. 0.) then + fteledifftrim = fmgytelediff(2,ix,1) else fteledifftrim = 0. endif # compute fmgyteleout if (isvacuummodel(1) .eq. 1) then - fmgyteleout(ix,1) = fmgytelemaxw(ix,1) + fmgyteleout(2,ix,1) = fmgytelemaxw(2,ix,1) elseif (isvacuummodel(1) .eq. 2) then - fmgyteleout(ix,1) = fteledifftrim + fmgyteleout(2,ix,1) = fteledifftrim elseif (isvacuummodel(1) .eq. 3) then - if (Kn(ix,1) <= Knb1) then - fmgyteleout(ix,1) = fteledifftrim - elseif (Kn(ix,1) >= Knb2) then - fmgyteleout(ix,1) = fmgytelemaxw(ix,1) + if (Kn(2,ix,1) <= Knb1) then + fmgyteleout(2,ix,1) = fteledifftrim + elseif (Kn(2,ix,1) >= Knb2) then + fmgyteleout(2,ix,1) = fmgytelemaxw(2,ix,1) else - fmgyteleout(ix,1) = (Kn(ix,1)-Knb1)/(Knb2-Knb1)*fmgytelemaxw(ix,1) + - . (Knb2-Kn(ix,1))/(Knb2-Knb1)*fteledifftrim + fmgyteleout(2,ix,1) = (Kn(2,ix,1)-Knb1)/(Knb2-Knb1)*fmgytelemaxw(2,ix,1) + + . (Knb2-Kn(2,ix,1))/(Knb2-Knb1)*fteledifftrim endif else - write (*,*) 'isvacuummodel > 3 not available, placeholder for future extention' + write (*,*) 'isvacuummodel > 3 not available for upg at outer wall, placeholder for future extention' endif # place holder for momentu influx fmgytelein - fmgytelein(ix,1) = 0. - yldot(iv2) = nurlxu * (fmiy(ix,ny,ifld) - (fmgyteleout(ix,1)-fmgytelein(ix,1)))/ + fmgytelein(2,ix,1) = 0. + yldot(iv2) = nurlxu * (fmiy(ix,ny,ifld) - (fmgyteleout(2,ix,1)-fmgytelein(2,ix,1)))/ . (vpnorm*sy(ix,ny)*fnorm(ifld)) else if (isupwoix(ix,ifld)==1) then #zero parallel momentum flux @@ -2881,37 +3018,38 @@ call sputchem (isch_sput(igsp), eng_sput, tvwallo(ix), if (istgonxy(ix,ny+1,igsp) == 1) then iv = idxtg(ix,ny+1,igsp) if (isvacuummodel(igsp) .gt. 0) then - fegytelemaxw(ix,igsp) = fngytelemaxw(ix,igsp)*2.*tg(ix,ny,igsp) - fegytelediff1(ix,igsp) = fngytelediff1(ix,igsp)*2.*tg(ix,ny,igsp) - fegytelediff2(ix,igsp) = fngytelediff2(ix,igsp)*2.5*tg(ix,ny,igsp) - fegytelediff3(ix,igsp) = fngytelediff3(ix,igsp)*2.5*tg(ix,ny,igsp) - fegytelediff(ix,igsp) = cfteleout*(fegytelediff1(ix,igsp) + fegytelediff2(ix,igsp) + fegytelediff3(ix,igsp)) + fegytelemaxw(2,ix,igsp) = fngytelemaxw(2,ix,igsp)*2.*tg(ix,ny,igsp) + fegytelediff1(2,ix,igsp) = fngytelediff1(2,ix,igsp)*2.*tg(ix,ny,igsp) + fegytelediff2(2,ix,igsp) = fngytelediff2(2,ix,igsp)*2.5*tg(ix,ny,igsp) + fegytelediff3(2,ix,igsp) = fngytelediff3(2,ix,igsp)*2.5*tg(ix,ny,igsp) + fegytelediff(2,ix,igsp) = cfteleout(2)*(fegytelediff1(2,ix,igsp) + + . fegytelediff2(2,ix,igsp) + fegytelediff3(2,ix,igsp)) # Trim diffusion to avoid negative influx - if (fngytelediff(ix,igsp) .gt. 0.) then - fteledifftrim = fegytelediff(ix,igsp) + if (fngytelediff(2,ix,igsp) .gt. 0.) then + fteledifftrim = fegytelediff(2,ix,igsp) else fteledifftrim = 0. endif # compute fegyteleout if (isvacuummodel(igsp) .eq. 1) then - fegyteleout(ix,igsp) = fegytelemaxw(ix,igsp) + fegyteleout(2,ix,igsp) = fegytelemaxw(2,ix,igsp) elseif (isvacuummodel(igsp) .eq. 2) then - fegyteleout(ix,igsp) = fteledifftrim + fegyteleout(2,ix,igsp) = fteledifftrim elseif (isvacuummodel(igsp) .eq. 3) then - if (Kn(ix,igsp) <= Knb1) then - fegyteleout(ix,igsp) = fteledifftrim - elseif (Kn(ix,igsp) >= Knb2) then - fegyteleout(ix,igsp) = fegytelemaxw(ix,igsp) + if (Kn(2,ix,igsp) <= Knb1) then + fegyteleout(2,ix,igsp) = fteledifftrim + elseif (Kn(2,ix,igsp) >= Knb2) then + fegyteleout(2,ix,igsp) = fegytelemaxw(2,ix,igsp) else - fegyteleout(ix,igsp) = (Kn(ix,igsp)-Knb1)/(Knb2-Knb1)*fegytelemaxw(ix,igsp) + - . (Knb2-Kn(ix,igsp))/(Knb2-Knb1)*fteledifftrim + fegyteleout(2,ix,igsp) = (Kn(2,ix,igsp)-Knb1)/(Knb2-Knb1)*fegytelemaxw(2,ix,igsp) + + . (Knb2-Kn(2,ix,igsp))/(Knb2-Knb1)*fteledifftrim endif else - write (*,*) 'isvacuummodel > 3 not available, placeholder for future extention' + write (*,*) 'isvacuummodel > 3 not available for tg at outer wall, placeholder for future extention' endif # place holder for momentu influx fmgytelein - fegytelein(ix,igsp) = 0. - yldot(iv) = nurlxg*( fegy(ix,ny,igsp) - (fegyteleout(ix,igsp)-fegytelein(ix,igsp)) )/ + fegytelein(2,ix,igsp) = 0. + yldot(iv) = nurlxg*( fegy(ix,ny,igsp) - (fegyteleout(2,ix,igsp)-fegytelein(2,ix,igsp)) )/ . (sy(ix,ny)*vpnorm*ennorm) else if (istgwcix(ix,igsp) == 0) then # fixed Tg