Skip to content

Commit

Permalink
Merge pull request #91 from LLNL/tg_rewrite_H2CRM-compatible_signfix
Browse files Browse the repository at this point in the history 
Fixes sign conventions and adds psorbgg
  • Loading branch information
@holm10
holm10 authored Nov 14, 2024
2 parents 9629d9f + d91c88b commit f31f2fd
Showing 1 changed file with 23 additions and 23 deletions.
46 changes: 23 additions & 23 deletions bbb/oderhs.m
View file
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

0 comments on commit f31f2fd

Please sign in to comment.