Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixes sign conventions and adds psorbgg #91

Merged
merged 2 commits into from
Nov 14, 2024
Merged

Fixes sign conventions and adds psorbgg #91

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f8b8ae6
Fixes sign conventions and adds psorbgg
holm10 Nov 14, 2024
d91c88b
Merge branch 'tg_rewrite_H2CRM-compatible' into tg_rewrite_H2CRM-comp…
holm10 Nov 14, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 23 additions & 23 deletions bbb/oderhs.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2350,18 +2350,20 @@ call xerrab('*** nhgsp must exceed 1 for ishymol=1 ***')
psorgc(ix,iy,2) = - ng(ix,iy,2)*nuiz(ix,iy,2)*vol(ix,iy) +
. psorbgg(ix,iy,2)
psorg(ix,iy,2) = psorgc(ix,iy,2) # no mol sor averaging
psordisg(ix,iy,2) = - ng(ix,iy,2)*nuiz(ix,iy,2)*vol(ix,iy)
psordis(ix,iy,2) = ng(ix,iy,2)*vol(ix,iy)*( 2*(1-ismolcrm)*
. ne(ix,iy)*(svdiss(te(ix,iy)) + sigvi_floor) +
. ismolcrm*cfcrma*sv_crumpet(te(ix,iy),ne(ix,iy),11))
psordisg(ix,iy,2) = ng(ix,iy,2)*nuiz(ix,iy,2)*vol(ix,iy)
. - psorbgg(ix,iy,2)
psordis(ix,iy,2) =
. - (1-ismolcrm)*2*psordisg(ix,iy,2)
. + ismolcrm*cfcrma*ng(ix,iy,2)*vol(ix,iy)
. * sv_crumpet(te(ix,iy),ne(ix,iy),11)
# 2 atoms per molecule in old model, rates from CRM for new
psordisg(ix,iy,1)=psordis(ix,iy,2)
psordis(ix,iy,1) = -cfcrmi*(2*psordisg(ix,iy,2)+
. psordis(ix,iy,2))
psor(ix,iy,1) = psor(ix,iy,1) + psordis(ix,iy,1)
c ... TODO: How to deal with diffusive atom model - is it maintained?
if(isupgon(1) .eq. 1) then
psor(ix,iy,iigsp) = psor(ix,iy,iigsp) + psordis(ix,iy,iigsp)
psor(ix,iy,iigsp) = psor(ix,iy,iigsp) - psordis(ix,iy,iigsp)
endif
enddo
enddo
Expand Down Expand Up @@ -4513,7 +4515,7 @@ c if (ifld .ne. iigsp) then
. eeli(ix,iy) = 13.6*ev +
. erliz(ix,iy)/(fac2sp*psor(ix,iy,1))

edisse(ix,iy)=-(1-ismolcrm)*ediss*ev*(0.5*psordis(ix,iy,2)) +
edisse(ix,iy)=-(1-ismolcrm)*ediss*ev*(-0.5*psordis(ix,iy,2)) +
. ismolcrm*ng(ix,iy,2)*vol(ix,iy)*
. sv_crumpet(te(ix,iy), ne(ix,iy), 20)
pradhyd(ix,iy)= ( (eeli(ix,iy)-ebind*ev)*psor(ix,iy,1)+
Expand Down Expand Up @@ -4625,7 +4627,7 @@ ccc call volave(nx, ny, j2, j5, i2, i5, ixp1(0,0), ixm1(0,0),
. + (1 + cftiexclg) * psicx(ix,iy) )

c Ion energy source from mol. diss ("Franck Condon")
seid(ix,iy) = cftiexclg * cfneut * cfneutsor_ei
seid(ix,iy) = - cftiexclg * cfneut * cfneutsor_ei
. * cnsor*eion*(1-ismolcrm)*ev*psordis(ix,iy,2)
. + emolia(ix,iy,1) + cftiexclg*emolia(ix,iy,2) # CRM FC

Expand Down Expand Up @@ -4653,7 +4655,7 @@ ccc call volave(nx, ny, j2, j5, i2, i5, ixp1(0,0), ixm1(0,0),

c Atom kinetic energy source from mol. diss
sead(ix,iy) = (1-ismolcrm)*(
. eion*ev*psordis(ix,iy,2)
. - eion*ev*psordis(ix,iy,2)
. )
. + emolia(ix,iy,2)

Expand Down Expand Up @@ -6049,7 +6051,7 @@ cc uu(ix,iy,iigsp) = uug(ix,iy,igsp)
. - fluxfacy*(fngy(ix,iy,igsp) - fngy(ix ,iy-1,igsp))
. + psgov_use(ix,iy,igsp)*vol(ix,iy)
if (igsp.eq.1 .and. ishymol.eq.1) resng(ix,iy,igsp) =
. resng(ix,iy,igsp)+psordis(ix,iy,2)
. resng(ix,iy,igsp)-psordis(ix,iy,2)
891 continue
892 continue
endif
Expand Down Expand Up @@ -6642,7 +6644,7 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
. + volpsorg(ix,iy,igsp)
. + psgov_use(ix,iy,igsp)*vol(ix,iy)
if (igsp.eq.1 .and. ishymol.eq.1)
. resng(ix,iy,igsp) = resng(ix,iy,igsp)+psordis(ix,iy,2)
. resng(ix,iy,igsp) = resng(ix,iy,igsp)-psordis(ix,iy,2)
resng(ix,iy,igsp) = resng(ix,iy,igsp) - cfneutdiv*
. cfneutdiv_fng*((fngx(ix,iy,igsp) - fngx(ix1,iy, igsp)) +
. fluxfacy*(fngy(ix,iy,igsp) - fngy(ix,iy-1,igsp)) )
Expand Down Expand Up @@ -7061,7 +7063,7 @@ cc uu(ix,iy,iigsp) = uug(ix,iy,igsp)
. + psgov_use(ix,iy,igsp)*vol(ix,iy)

if (igsp.eq.1 .and. ishymol.eq.1)
. resng(ix,iy,igsp) = resng(ix,iy,igsp)+psordis(ix,iy,2)
. resng(ix,iy,igsp) = resng(ix,iy,igsp)-psordis(ix,iy,2)
891 continue
892 continue
endif
Expand Down Expand Up @@ -7479,7 +7481,7 @@ call fd2tra (nx,ny,floxg,floyg,conxg,conyg,
. - fngy(ix,iy,igsp) + fngy(ix ,iy-1,igsp)
. + psgov_use(ix,iy,igsp)*vol(ix,iy)
if (igsp.eq.1 .and. ishymol.eq.1)
. resng(ix,iy,igsp) = resng(ix,iy,igsp)+psordis(ix,iy,2)
. resng(ix,iy,igsp) = resng(ix,iy,igsp)-psordis(ix,iy,2)
891 continue
892 continue
else if (isupgon(igsp).eq.1) then
Expand Down Expand Up @@ -7880,9 +7882,7 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp),

* Internal molecular energy sink due to dissociation
* ----------------------------------------------------
reseg(ix,iy,2) = reseg(ix,iy,2)
. + (1-ismolcrm) *psorg(ix,iy,2)*1.5*tg(ix,iy,2)
. + ismolcrm * 1.5*tg(ix,iy,2)*psordisg(ix,iy,2)
reseg(ix,iy,2) = reseg(ix,iy,2) + psorg(ix,iy,2)*1.5*tg(ix,iy,2)


* Drift heating energy source for molecules
Expand All @@ -7906,34 +7906,34 @@ c are corrections for (v_m - v_a)**2. However, the original
v2cc = (cfnidhg2**0.5)*0.5*(v2(ix,iy,iigsp)+v2(ix1,iy,iigsp))

reseg(ix,iy,1) = reseg(ix,iy,1)
. + cfnidh*cfnidhdis*0.5*mg(1)* (upgcc**2 + vycc**2 + v2cc**2)
. - cfnidh*cfnidhdis*0.5*mg(1)* (upgcc**2 + vycc**2 + v2cc**2)
. * psordis(ix,iy,2)

seic(ix,iy) = seic(ix,iy)
. + cftiexclg * cfneut * cfneutsor_ei * cnsor * cfnidhdis
. - cftiexclg * cfneut * cfneutsor_ei * cnsor * cfnidhdis
. * 0.5*mg(1)*(upgcc**2 + vycc**2 + v2cc**2)
. * ( (1-ismolcrm)*psordis(ix,iy,2) + ismolcrm*psordis(ix,iy,1) )

uuxgcc = (cfnidhmol**0.5)*0.5*(uuxg(ix,iy,2)+uuxg(ix1,iy,2))**2
vygcc = (cfnidhmol**0.5)*0.5*(vyg(ix,iy,2)+vyg(ix1,iy,2))**2
v2gcc = 0. #.. molecule v in the tol direction, it seems to be assumed as 0 in neudifpg?
reseg(ix,iy,1) = reseg(ix,iy,1)
. + cfnidhdis*0.5*mg(1)*(uuxgcc**2 + vygcc**2 + v2gcc**2 )*psordis(ix,iy,2)
. - cfnidhdis*0.5*mg(1)*(uuxgcc**2 + vygcc**2 + v2gcc**2 )*psordis(ix,iy,2)

seic(ix,iy) = seic(ix,iy)
. + cftiexclg*cfnidhdis*0.5*mg(1)*(uuxgcc**2 + vygcc**2 + v2gcc**2 )
. * ( (1-ismolcrm)*psordis(ix,iy,2) + ismolcrm*psordis(ix,iy,1) )
. - ( (1-ismolcrm)*psordis(ix,iy,2) + ismolcrm*psordis(ix,iy,1) )

uuxgcc = cfnidhmol*0.25*(uuxg(ix,iy,2)+uuxg(ix1,iy,2))
. *(uuxg(ix,iy,1)+uuxg(ix1,iy,1))
vygcc = cfnidhmol*0.25*(vyg(ix,iy,2)+vyg(ix1,iy,2))
. *(vyg(ix,iy,1)+vyg(ix1,iy,1))
v2gcc = 0.
reseg(ix,iy,1) = reseg(ix,iy,1)
. - cfnidhdis*mg(1)*(uuxgcc + vygcc + v2gcc)*psordis(ix,iy,2)
. + cfnidhdis*mg(1)*(uuxgcc + vygcc + v2gcc)*psordis(ix,iy,2)

seic(ix,iy) = seic(ix,iy)
. - cftiexclg*cfnidhdis*mg(1)*(uuxgcc + vygcc + v2gcc)
. + cftiexclg*cfnidhdis*mg(1)*(uuxgcc + vygcc + v2gcc)
. * ( (1-ismolcrm)*psordis(ix,iy,2) + ismolcrm*psordis(ix,iy,1) )
* Start new Molecular Drift heating implementation
* ----------------------------------------------------
Expand All @@ -7947,12 +7947,12 @@ c are corrections for (v_m - v_a)**2. However, the original
v2gcc = 0. #.. molecule v in the tol direction, it seems to be assumed as 0 in neudifpg?

reseg(ix,iy,1) = reseg(ix,iy,1)
. + cfnidhdis*0.5*mg(1)
. - cfnidhdis*0.5*mg(1)
. * ((uuxgcc-upgcc)**2 + (vygcc-vycc)**2 + (v2gcc-v2cc)**2 )
. * psordis(ix,iy,2)

seic(ix,iy) = seic(ix,iy)
. + cftiexclg*cfnidhdis*0.5*mg(1)
. - cftiexclg*cfnidhdis*0.5*mg(1)
. * ((uuxgcc-upgcc)**2 + (vygcc-vycc)**2 + (v2gcc-v2cc)**2 )
. * ( (1-ismolcrm)*psordis(ix,iy,2) + ismolcrm*psordis(ix,iy,1) )
endif # End new drift heating implementation
Expand Down
Loading