Skip to content

Commit

Permalink
Fixes Ta equation when molecules included
Browse files Browse the repository at this point in the history
  • Loading branch information
@holm10
holm10 committed Nov 16, 2024
1 parent 4e28575 commit 05350b2
Showing 1 changed file with 25 additions and 27 deletions.
52 changes: 25 additions & 27 deletions bbb/oderhs.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4658,7 +4658,7 @@ ccc call volave(nx, ny, j2, j5, i2, i5, ixp1(0,0), ixm1(0,0),
. * (psorrg(ix,iy,1)+psicx(ix,iy))

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

Expand Down Expand Up @@ -7848,12 +7848,12 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp),

* -- total energy residual and equipartition --

do iy = j2, j5
iy1 = max(0,iy-1)
do ix = i2, i5
ix1 = ixm1(ix,iy)
do igsp = 1, ngsp
do iy = j2, j5
iy1 = max(0,iy-1)
do ix = i2, i5
ix1 = ixm1(ix,iy)

do igsp = 1, ngsp
* Compute thermal equipartition rate with ion-atom fluid
* ------------------------------------------------------------
if (nisp >= 2) then # uses ni(,,2), so must have atoms
Expand All @@ -7866,13 +7866,14 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp),
* GAS-ION/ATOM TERMS
* ------------------------------------------------------------


* Divergence of gaseous flows & v-grad-P heating
* ------------------------------------------------------------
reseg(ix,iy,igsp)= -( fegx(ix,iy,igsp)-fegx(ix1,iy, igsp)+
. fegy(ix,iy,igsp)-fegy(ix, iy1,igsp) )
. + segc(ix,iy,igsp)



* ------------------------------------------------------------
* GAS-GAS TERMS
* ------------------------------------------------------------
Expand Down Expand Up @@ -7902,14 +7903,9 @@ call fd2tra (nx,ny,floxge(0:nx+1,0:ny+1,igsp),
reseg(ix,iy,1) = reseg(ix,iy,1)
. - cftgeqp*vol(ix,iy)*kelighg(igsp)
. * ng(ix,iy,igsp)*ng(ix,iy,1)*(tg(ix,iy,1)-tg(ix,iy,igsp))
endif


* ----------------------------------------------------
* COUPLE MOLECULES TO ATOMS AND IONS
* ----------------------------------------------------
if (ishymol .eq. 1 .and. igsp .eq. 2) then
c Add check for inertial atoms?
if (ishymol.eq.1 .and. igsp.eq.2) then #..D2 dissociation

* Internal ion/atom energy source due to dissociation
* ----------------------------------------------------
Expand Down Expand Up @@ -7948,26 +7944,26 @@ c are corrections for (v_m - v_a)**2. However, the original

c Only apply drift heating for inertial atoms?
* ----------------------------------------------------
if (1 .eq. 1) then
if (1.eq.1) then

upgcc = 0.5*(up(ix,iy,iigsp)+up(ix1,iy,iigsp))
vycc = (cfnidhgy**0.5)*0.5*(vy(ix,iy,iigsp)+vy(ix1,iy,iigsp))
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)
. * psordis(ix,iy,2)

seic(ix,iy) = seic(ix,iy)
. - 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
uuxgcc = (cfnidhmol**0.5)*0.5*(uuxg(ix,iy,2)+uuxg(ix1,iy,2))
vygcc = (cfnidhmol**0.5)*0.5*(vyg(ix,iy,2)+vyg(ix1,iy,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)

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) )
Expand All @@ -7977,21 +7973,20 @@ c are corrections for (v_m - v_a)**2. However, the original
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)

reseg(ix,iy,1) = reseg(ix,iy,1)
. + 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
* ----------------------------------------------------

else
upgcc = 0.5*(up(ix,iy,iigsp)+up(ix1,iy,iigsp))
vycc = (cfnidhgy**0.5)*0.5*(vy(ix,iy,iigsp)+vy(ix1,iy,iigsp))
v2cc = (cfnidhg2**0.5)*0.5*(v2(ix,iy,iigsp)+v2(ix1,iy,iigsp))

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
uuxgcc = (cfnidhmol**0.5)*0.5*(uuxg(ix,iy,2)+uuxg(ix1,iy,2))
vygcc = (cfnidhmol**0.5)*0.5*(vyg(ix,iy,2)+vyg(ix1,iy,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)
Expand All @@ -8003,16 +7998,19 @@ c are corrections for (v_m - v_a)**2. However, the original
. - 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
endif
endif


endif
endif
#..zml place holder for neutral-neutral collision,
#.. not included above?
enddo
enddo
enddo



if((isudsym==1.or.geometry.eq.'dnXtarget') .and. nxc > 1) then
do igsp = 1, ngsp
do iy = j2, j5
Expand Down

0 comments on commit 05350b2

Please sign in to comment.