Combine H2 CRM branch with Tg branch

.github/workflows/build-test-wheel.yml

@@ -30,16 +30,17 @@ jobs:
     - name: Return to pytests branch
       uses: actions/checkout@v4
       with:
-        ref: pytests
+        ref: pytests-develop
         repository: LLNL/UEDGE 
     - name: Test with pytest
       run: |
         # Create a clean UEDGE instance for every test
-        pytest --isolate --tb=native pytests/${{ matrix.tests }}
+        pytest --isolate --tb=native -s pytests/${{ matrix.tests }}
 
 
   build_wheels:
     name: Build wheels on ${{ matrix.os }} for Python ${{ matrix.pyver }}
+    if: false
  #   if: github.event_name == 'pull_request'
     needs: run-all-tests
     runs-on: ${{ matrix.os }}
@@ -86,6 +87,7 @@ jobs:
 
   build_wheels_pypy:
     name: Build PyPy wheels on ${{ matrix.os }} 
+    if: false
 #    if: github.event_name == 'pull_request'
     needs: run-all-tests
     runs-on: ${{ matrix.os }}
@@ -193,6 +195,7 @@ jobs:
 
   gather_wheels:
     name: Gather built wheels
+    if: false
     needs: [build_wheels, build_wheels_pypy]
 #    if: github.event_name == 'pull_request'
     runs-on: ubuntu-latest

aph/aph.v

@@ -13,6 +13,8 @@ sgvcxc   real [m^2/s] /2.e-14/# const. value of sigv_cx for issgvcxc=1;
 isaphdir integer       /1/    #=1 uses aphdir; =0 uses explicit rate file names
 crmnfname character*120 /"crumpet_nrates.dat"/ # Name of CRM rate data file
 crmefname character*120 /"crumpet_Erates.dat"/ # Name of CRM E-sink data file
+crmdir character*120 /"."/  # Path to folder containing CRM data files: data files
+                            # specified by crmnfname and crmename
 aphdir	character*120 # name of directory containing data files
 data_directory	character*120 # another dirname containing data files. This is to be be passed in
 

aph/aphrates.m

@@ -1379,7 +1379,7 @@ real function sv_crumpet (te, ne, rate)
             sv_crumpet=0
       else
 c     compute abscissae --
-         zloge=log(te/ev)
+         zloge=log(te/ev_aph)
          rle=max(crlemin, min(zloge,crlemax))
          zlogd=log10(ne)
          rld=max(crldmin, min(zlogd,crldmax))
@@ -1461,7 +1461,7 @@ call remark("Rate option not recognized!")
                 jd=jd-1
                 je=je-1
                 fjd=((ne/1e6)-(10**(10+0.5*jd)))/ (10**(10+0.5*jd)*(10**0.5 - 1) )
-                fje=(te/ev - (10**(-1.2+0.1*je)))/( (10**(-1.2+0.1*je))*(10**0.1-1))
+                fje=(te/ev_aph - (10**(-1.2+0.1*je)))/( (10**(-1.2+0.1*je))*(10**0.1-1))
 
                 svd_crm1=svd_crm11+fjd*(svd_crm12-svd_crm11)
                 svd_crm2=svd_crm21+fjd*(svd_crm22-svd_crm21)

aph/aphread.m

@@ -136,15 +136,17 @@ call readehr2(TRIM(adname))
       subroutine crumpetread
 
 c ... Set up tables for hydrogenic molecular processes using a CRUMPET
-c ... CRM file
+c ... CRM file. Looks for files crmnfname and crmefname in the directory
+c ... specified by crmdir. If not found, looks for the file in aphdir,
+c ... and finally in the CWD
 
       implicit none
 
 c ... Common blocks:
       Use(Dim)
       Use(UEpar)        # ismolcrm
       Use(Share)        # nhdf, hdfilename
-      Use(Data_input)   # aphdir, data_directory, crmefname, crmnfname
+      Use(Data_input)   # aphdir, data_directory, crmdir, crmefname, crmnfname
       Use(Rtcrumpet)      # cmpd, cmpe
 
 c ... Function:
@@ -187,10 +189,11 @@ call gallot("Rtcrumpet",0)
          enddo
 
                 if (ismolcrm .ne. 0) then
-
-                     call findFile(crmnfname, aphdirx, dataDir, adname, isaphdir)
+c                     call findFile(crmnfname, aphdirx, aphdirx, adname, isaphdir)
+                     call findFile(crmnfname, crmdir, aphdirx, adname, isaphdir)
                      call readcrumpetn(TRIM(adname))
-                     call findFile(crmefname, aphdirx, dataDir, adname, isaphdir)
+                     call findFile(crmefname, aphdirx, aphdirx, adname, isaphdir)
+c                     call findFile(crmefname, crmdir, aphdirx, adname, isaphdir)
                      call readcrumpete(TRIM(adname))
 
 
@@ -969,6 +972,7 @@ subroutine readcrumpetn (fname)
 c     local variables --
       integer ios, nget, jd, jt
       character*80 zdummy
+      character(len=500) :: io_emsg
       real crmdummy(cmpe,cmpd)
 
 c     procedures --
@@ -979,12 +983,11 @@ real crmdummy(cmpe,cmpd)
 c----------------------------------------------------------------------c
 
       call freeus(nget)
-      open (nget, file=fname, form='formatted', iostat=ios,
-     .     status='old')
-      if (ios .ne. 0) then
-      write(*,*) fname
-         call xerrab('**** CRM rate file not found; set aphdir
-     . path')
+      open (nget, file=TRIM(fname), form='formatted', iostat=ios,
+     .     status='old',iomsg=io_emsg)
+      if (ios .ne. 0) then 
+         print *, trim(io_emsg)
+         call xerrab("")
       endif
 
 c     Atom depletion rate in CRM (cm**3/sec):
@@ -1053,6 +1056,7 @@ subroutine readcrumpete (fname)
 c     local variables --
       integer ios, nget, jd, jt
       character*80 zdummy
+      character(len=500) :: io_emsg
       real crmdummy(cmpe,cmpd)
 
 c     procedures --
@@ -1064,10 +1068,10 @@ real crmdummy(cmpe,cmpd)
 
       call freeus(nget)
       open (nget, file=fname, form='formatted', iostat=ios,
-     .     status='old')
+     .     status='old',iomsg=io_emsg)
       if (ios .ne. 0) then
-         call xerrab('**** CRM rate file not found; set aphdir
-     . path')
+         print *, trim(io_emsg)
+         call xerrab("")
       endif
 
 c     Electron energy change due to atoms (W cm**3):
@@ -1180,11 +1184,11 @@ c     Atom radiation source due to external interactions (W cm**3):
 c     convert to SI units:
       do jt=1,cmpe
          do jd=1,cmpd
-            crmselm(jt,jd)= crmselm(jt,jd)*ev
-            crmsiam(jt,jd)= crmsiam(jt,jd)*ev
-            crmspotm(jt,jd)= crmspotm(jt,jd)*ev
-            crmsrada(jt,jd)= crmsrada(jt,jd)*ev
-            crmsradm(jt,jd)= crmsradm(jt,jd)*ev
+            crmselm(jt,jd)= crmselm(jt,jd)*ev_aph
+            crmsiam(jt,jd)= crmsiam(jt,jd)*ev_aph
+            crmspotm(jt,jd)= crmspotm(jt,jd)*ev_aph
+            crmsrada(jt,jd)= crmsrada(jt,jd)*ev_aph
+            crmsradm(jt,jd)= crmsradm(jt,jd)*ev_aph
          enddo
       enddo
 

bbb/bbb.v

@@ -2223,8 +2223,8 @@ psorbgg(0:nx+1,0:ny+1,1:ngsp) _real  [part/s]  # diag artific neut backg source
 psorbgz(0:nx+1,0:ny+1)        _real  [part/s]  # diag artific impur backg source
 erliz(0:nx+1,0:ny+1)          _real  [J/s]     # H rad'n loss for ioniz'n
 edisse(0:nx+1,0:ny+1)         _real  [J/s]     # Elec E loss due to mol interactions
-emolia(0:nx+1,0:ny+1)         _real  [J/s]     # i/a E change due to mol interactions
-eiamoldiss(0:nx+1,0:ny+1)     _real  [J/s]     # i/a enegy density incr, mol diss
+emolia(0:nx+1,0:ny+1,1:nisp)  _real  [J/s]     # i/a E change due to mol interactions
+eiamoldiss(0:nx+1,0:ny+1,1:nisp)    _real  [J/s]     # i/a enegy density incr, mol diss
 erlrc(0:nx+1,0:ny+1)          _real  [J/s]     # H rad'n loss for recom'n
 vsoreec(0:nx+1,0:ny+1)	      _real  [J/s]     # cell ctr tot elec vol eng source
 vsoree(0:nx+1,0:ny+1)	      _real  [J/s]     # cell ave tot elec vol eng source
@@ -2242,6 +2242,20 @@ seec(0:nx+1,0:ny+1)           _real
 seev(0:nx+1,0:ny+1)           _real
 seic(0:nx+1,0:ny+1)           _real
 seiv(0:nx+1,0:ny+1)           _real
+seik(0:nx+1,0:ny+1)           _real +work   # Kinetic energy source from recom.
+                                            # and ioniz. (ions)
+seid(0:nx+1,0:ny+1)           _real +work   # Kinetic energy source from 
+                                            # dissociation (ions)
+seidh(0:nx+1,0:ny+1)          _real +work   # Drift heating (ions)
+seit(0:nx+1,0:ny+1)           _real +work   # Internal energy source/sink
+                                            # from ioniz and recom (ions)
+psicx(0:nx+1,0:ny+1)          _real +work   # CX rate (ions)
+seak(0:nx+1,0:ny+1)           _real +work   # Kinetic energy sink/source from 
+                                            # rec and CX (atoms)
+sead(0:nx+1,0:ny+1)           _real +work   # Kinetic energy source from 
+                                            # dissociation (atoms)
+seadh(0:nx+1,0:ny+1)          _real +work   # Drift heating (atoms)
+
 segc(0:nx+1,0:ny+1,1:ngsp)    _real [J/(sm**3)]#v_grad_P for neutral eng. eqn
 resco(0:nx+1,0:ny+1,1:nisp)   _real
 resng(0:nx+1,0:ny+1,1:ngsp)   _real

bbb/boundary.m

@@ -18,7 +18,7 @@ real yl(neq), yldot(neq)
       Use(Phyvar)
       Use(UEpar)    # isnewpot,r0slab,cslim,dcslim,csfaclb,csfacrb,csfacti,
                     # isnion,isupon,isteon,istion,isngon,isnionxy,isuponxy,
-                    # isteonxy,istionxy,isngonxy,isphionxy
+                    # isteonxy,istionxy,isngonxy,isphionxy, ismolcrm
       Use(Aux)      # ixmp
       Use(Coefeq)   # fac2sp,cf2ef,exjbdry
       Use(Bcond)    # iflux,ncore,tcoree,tcorei,tbmin,nbmin,ngbmin,
@@ -1731,7 +1731,7 @@ ccc              yldot(iv2) = nurlxp*(phiwo(ix) - phi(ix,ny))/temp0
                   if (ishymol.eq.1 .and. igsp.eq.2) then
                     ta0 = engbsr * max(tg(1,iy,1),temin*ev)
                     vxa = 0.25 * sqrt( 8*ta0/(pi*mg(1)) )
-                    flxa = (1-alblb(iy,1,1))*ng(1,iy,1)*vxa*sx(0,iy)
+                    flxa = ismolcrm*(1-alblb(iy,1,1))*ng(1,iy,1)*vxa*sx(0,iy)
 
                     if (isupgon(1) .eq. 1) then  # two atoms per molecule
                       flux_inc = 0.5*( fnix(0,iy,1) + fnix(0,iy,2) + flxa)
@@ -2086,7 +2086,7 @@ cc               if (recyce .le. 0) bcen = 0.  # gets back to old case
                if (ishymol.eq.1 .and. igsp.eq.2) then
                  ta0 = max(tg(ixt1,iy,1), temin*ev)
                  vxa = 0.25 * sqrt( 8*ta0/(pi*mg(1)) )
-                 flxa = (1-alblb(iy,1,jx))*ng(ixt1,iy,1)*vxa*sx(ixt,iy)
+                 flxa = ismolcrm*(1-alblb(iy,1,jx))*ng(ixt1,iy,1)*vxa*sx(ixt,iy)
 
                  if (isupgon(1) .eq. 1) then  # two atoms for one molecule
                    flux_inc = 0.5*( fnix(ixt,iy,1) + fnix(ixt,iy,2) +flxa) 
@@ -2405,7 +2405,7 @@ ccc          yldot(iv) = -nurlxp*(phi(0,iy)-phi(1,iy))/temp0
                   if (ishymol.eq.1 .and. igsp.eq.2) then
                     ta0 = engbsr * max(tg(nx,iy,1),temin*ev)
                     vxa = 0.25 * sqrt( 8*ta0/(pi*mg(1)) )
-                    flxa=(1-albrb(iy,1,nxpt))*ng(nx,iy,1)*vxa*sx(nx,iy)
+                    flxa= ismolcrm*(1-albrb(iy,1,nxpt))*ng(nx,iy,1)*vxa*sx(nx,iy)
                     if (isupgon(1) .eq. 1) then  # two atoms for one molecule
                       flux_inc = 0.5*( fnix(nx,iy,1) + fnix(nx,iy,2) -flxa) 
                     else
@@ -2777,7 +2777,7 @@ cc               if (recyce .le. 0) bcen = 0.  # gets back to old case
                if (ishymol.eq.1 .and. igsp.eq.2) then
                 ta0 = max(tg(ixt1,iy,1), temin*ev)
                 vxa = 0.25 * sqrt( 8*ta0/(pi*mg(1)) )
-                flxa=(1-albrb(iy,1,jx))*ng(ixt1,iy,1)*vxa*sx(ixt1,iy)
+                flxa= ismolcrm*(1-albrb(iy,1,jx))*ng(ixt1,iy,1)*vxa*sx(ixt1,iy)
 
                  if (isupgon(1) .eq. 1) then  # two atoms for one molecule
                    flux_inc = 0.5*( fnix(ixt1,iy,1) +fnix(ixt1,iy,2)-flxa)