From 1418647cb2ab9fb836c13152df67c14a92e5e0f1 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Thu, 21 Nov 2024 10:40:30 -0800 Subject: [PATCH] Replaces radial maxwellian fluxes w/ function calls --- bbb/boundary.m | 41 ++++++++++++++++++++--------------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/bbb/boundary.m b/bbb/boundary.m index 757c2832..d061601f 100755 --- a/bbb/boundary.m +++ b/bbb/boundary.m @@ -244,30 +244,29 @@ call xerrab("*** isngcore=2 option unvailable ***") endif else # ix is not part of the core boundary if (iscpli(ix) .eq. 1) call wsmodi(1) - t0 = max(tg(ix,0,1),tgmin*ev) - vyn = sqrt( 0.5*t0/(pi*mi(1)) ) fng_chem = 0. - do ii = 1, ngsp #chem sputt of hydrogen - strange = 0 - flx_incid = ng(ix,1,ii)*0.25*sqrt(8*t0/(pi*mg(ii))) - fng_chem= fng_chem + chemsputi(1,ii)*flx_incid* - . sy(ix,0) + do ii = 1, ngsp #chem sputt of hydrogen - strange = 0 + fng_chem= fng_chem + chemsputi(1,ii)*onesided_maxwellian( + . tg(ix,0,1), ng(ix,1,ii), mi(1), sy(ix,0), tgmin*ev) enddo - 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) - yldot(iv1) = -nurlxg*( fniy(ix,0,ifld) + fng_alb - - . fng_chem ) / - . (vyn*sy(ix,0)* n0(ifld)) -cc if (fng_chem .ne. 0.) yldot(iv1) = -nurlxg*( -cc . fniy(ix,0,ifld) - fng_chem )/(vyn*sy(ix,0)*n0(ifld)) + yldot(iv1) = -nurlxg*( fniy(ix,0,ifld) + (1-albedoi(ix,1))* + . onesided_maxwellian( + . tg(ix,0,1), harmave(ni(ix,0,1), ni(ix,1,ifld)), + . mi(1), sy(ix,0), tgmin*ev) -fng_chem ) + . / onesided_maxwellian( tg(ix,0,1), n0(ifld), + . mi(1), sy(ix,0), tgmin*ev) c... Caution: the wall source models assume gas species 1 only is inertial if(matwalli(ix) .gt. 0) then if (recycwit(ix,1,1) .gt. 0.) then fniy_recy = recycwit(ix,1,1)*fac2sp*fniy(ix,0,1) if (isrefluxclip==1) fniy_recy=min(fniy_recy,0.) yldot(iv1)=-nurlxg*( fniy(ix,0,ifld) + fniy_recy - - . fngyi_use(ix,1) - fngysi(ix,1) + fng_alb - - . fng_chem ) / (vyn*n0(ifld)*sy(ix,0)) + . fngyi_use(ix,1) - fngysi(ix,1) + (1-albedoi(ix,1)) + . * onesided_maxwellian( + . tg(ix,0,1), harmave(ni(ix,0,ifld), ni(ix,1,ifld)), + . mi(1), sy(ix,0), tgmin*ev) - fng_chem ) + . / onesided_maxwellian( tg(ix,0,1), n0(ifld), + . mi(1), sy(ix,0), tgmin*ev) elseif (recycwit(ix,1,1) < -1) then yldot(iv1)=nurlxg*(ngbackg(1)-ni(ix,0,ifld))/n0(ifld) elseif (recycwit(ix,1,1) .le. 0.) then # treat recycwit as albedo @@ -303,7 +302,7 @@ cc if (fng_chem .ne. 0.) yldot(iv1) = -nurlxg*( yldot(iv1) = nurlxn*(ncore(ifld)-ni(ix,0,ifld))/ . n0(ifld) endif - elseif (isnicore(ifld) .eq. 3) then # const ni; flux set to + elseif (isnicore(ifld) .eq. 3) then # const ni; flux set to # curcore-recycc*fngy yldot(iv1) = -nurlxn*( ni(ix,0,ifld) - . ni(ixp1(ix,0),0,ifld) ) / n0(ifld) @@ -318,7 +317,7 @@ cc if (fng_chem .ne. 0.) yldot(iv1) = -nurlxg*( fngytotc = fngytotc + fngy(ii,0,1) #generalize if (ii==ix_fl_bc) break enddo - if (ifld > 1) fngytotc = 0. # crude fix for imp. + if (ifld > 1) fngytotc = 0. # crude fix for imp. yldot(iv1)= -nurlxn*(qe*fniytotc + . qe*recycc(1)*fngytotc - . curcore(ifld))/(qe*vpnorm*n0(ifld)) @@ -335,7 +334,7 @@ cc if (fng_chem .ne. 0.) yldot(iv1) = -nurlxg*( . ncoremin(ifld) ) / n0(ifld) endif else - call xerrab ('** isnicore value not valid option **') + call xerrab ('** isnicore value not valid option **') endif elseif (isnwconiix(ix,ifld) .eq. 0) then # ix not on core boundary; set zero gradient (or flux) @@ -371,13 +370,13 @@ call xerrab ('** isnicore value not valid option **') do jx = 1, nxpt if(isfixlb(jx).ne.2 .and. ixmnbcl*iymnbcl.eq.1 .and. . isnionxy(ixlb(jx),0,ifld)==1) then - yldot(idxn(ixlb(jx),0,ifld)) = nurlxn * + yldot(idxn(ixlb(jx),0,ifld)) = nurlxn * . ( ave(ni(ixlb(jx),1,ifld),ni(ixlb(jx)+1,0,ifld)) . - ni(ixlb(jx),0,ifld) ) / n0(ifld) endif if (isfixrb(jx).ne.2 .and. ixmxbcl*iymnbcl.eq.1 .and. . isnionxy(ixrb(jx)+1,0,ifld)==1) then - yldot(idxn(ixrb(jx)+1,0,ifld)) = nurlxn * + yldot(idxn(ixrb(jx)+1,0,ifld)) = nurlxn * . ( ave(ni(ixrb(jx)+1,1,ifld),ni(ixrb(jx),0,ifld)) . - ni(ixrb(jx)+1,0,ifld) ) / n0(ifld) endif