Merge branch 'development'

Author: zingale

.zenodo.json

@@ -21,6 +21,10 @@
             "name": "Katz, Maximilian",
             "orcid": "0000-0003-0439-4556"
         },
+        {   "affiliation": "Department of Physics and Astronomy, Stony Brook University",
+            "name": "Li, Xinlong",
+            "orcid": "0000-0002-3143-6136"
+        },
         {
             "affiliation": "Los Alamos National Laboratory",
             "name": "Malone, Chris",

CHANGES

@@ -1,4 +1,10 @@
-# 19.02
+# 19.05
+
+  -- The aprox21 network was missing the analytic Jacobian term for
+     the derivative of He4 with respect to Ni56. This is fixed. (#175)
+
+  -- The numerical Jacobian module, used by the BS and VBDF integrators
+     had some wrong scalings.  These have now been fixed (#179, #180)
 
 # 19.01
 

EOS/helmholtz/actual_eos.F90

@@ -60,7 +60,7 @@ module actual_eos_module
     integer, parameter          :: max_newton = 100
 
     ! 2006 CODATA physical constants
-private
+
     ! Math constants
     double precision, parameter :: pi       = 3.1415926535897932384d0
 

integration/VODE90/cuVODE/linear_algebra_modules/blas_module.F90

@@ -4,62 +4,6 @@ module blas_module
 
 contains
 
-  SUBROUTINE DCOPYN(N,DX,INCX,DY,INCY)
-    !$gpu
-  ! Only operates on arrays of size N
-
-    INTEGER INCX,INCY,N
-    DOUBLE PRECISION DX(N),DY(N)
-! *  Purpose
-! *  =======
-! *
-! *     copies a array, x, to a array, y.
-! *     uses unrolled loops for increments equal to one.
-! *     jack dongarra, linpack, 3/11/78.
-! *     modified 12/3/93, array(1) declarations changed to array(*)
-    INTEGER I,IX,IY,M,MP1
-
-    IF (N.LE.0) RETURN
-    IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
-
-    ! code for unequal increments or equal increments
-    ! not equal to 1
-
-    IX = 1
-    IY = 1
-    IF (INCX.LT.0) IX = (-N+1)*INCX + 1
-    IF (INCY.LT.0) IY = (-N+1)*INCY + 1
-    DO I = 1,N
-       DY(IY) = DX(IX)
-       IX = IX + INCX
-       IY = IY + INCY
-    end do
-    RETURN
-
-    ! code for both increments equal to 1
-
-    ! clean-up loop
-
-20  M = MOD(N,7)
-    IF (M.EQ.0) GO TO 40
-    DO I = 1,M
-       DY(I) = DX(I)
-    end do
-    IF (N.LT.7) RETURN
-40  MP1 = M + 1
-    DO I = MP1,N,7
-       DY(I) = DX(I)
-       DY(I+1) = DX(I+1)
-       DY(I+2) = DX(I+2)
-       DY(I+3) = DX(I+3)
-       DY(I+4) = DX(I+4)
-       DY(I+5) = DX(I+5)
-       DY(I+6) = DX(I+6)
-    end do
-    RETURN
-  end SUBROUTINE DCOPYN
-
-  
   SUBROUTINE DCOPY(N,DX,INCX,DY,INCY)
     !$gpu
     INTEGER INCX,INCY,N

integration/VODE90/cuVODE/source/cuvode.F90

@@ -302,7 +302,7 @@ subroutine dvode(vstate, rwork, IWORK, ITASK, IOPT, MF)
     vstate % JSTART = -1
     IF (vstate % NQ .LE. VODE_MAXORD) GO TO 90
     ! MAXORD was reduced below NQ.  Copy YH(*,MAXORD+2) into SAVF. ---------
-    CALL DCOPYN(VODE_NEQS, rwork % wm, 1, rwork % savf, 1)
+    rwork % savf(1:VODE_NEQS) = rwork % wm(1:VODE_NEQS)
 
     ! Reload WM(1) = RWORK % wm(1), since LWM may have changed. ---------------
 90  continue
@@ -356,7 +356,7 @@ subroutine dvode(vstate, rwork, IWORK, ITASK, IOPT, MF)
     CALL f_rhs (vstate % T, vstate % Y, rwork % yh(:,2), vstate % RPAR)
     vstate % NFE = 1
     ! Load the initial value array in YH. ---------------------------------
-    CALL DCOPYN(VODE_NEQS, vstate % Y, 1, rwork % YH(:,1), 1)
+    rwork % YH(1:VODE_NEQS,1) = vstate % Y(1:VODE_NEQS)
 
     ! Load and invert the EWT array.  (H is temporarily set to 1.0.) -------
     vstate % NQ = 1
@@ -641,7 +641,7 @@ subroutine dvode(vstate, rwork, IWORK, ITASK, IOPT, MF)
     ! -----------------------------------------------------------------------
 
 400 CONTINUE
-    CALL DCOPYN(VODE_NEQS, rwork % YH(:,1), 1, vstate % Y, 1)
+    vstate % Y(1:VODE_NEQS) = rwork % YH(1:VODE_NEQS,1)
 
     vstate % T = vstate % TN
     IF (ITASK .NE. 4 .AND. ITASK .NE. 5) GO TO 420
@@ -732,7 +732,7 @@ subroutine dvode(vstate, rwork, IWORK, ITASK, IOPT, MF)
     IWORK(16) = IMXER
     ! Set Y array, T, and optional output. --------------------------------
 580 CONTINUE
-    CALL DCOPYN(VODE_NEQS, rwork % YH(:,1), 1, vstate % Y, 1)
+    vstate % Y(1:VODE_NEQS) = rwork % YH(1:VODE_NEQS,1)
 
     vstate % T = vstate % TN
     IWORK(11) = vstate % NST

integration/VODE90/cuVODE/source/cuvode_dacopy.F90

@@ -39,7 +39,7 @@ subroutine dacopy(NROW, NCOL, A, NROWA, B, NROWB)
     !$gpu
 
     do IC = 1,NCOL
-       CALL DCOPYN (NROW, A(:,IC), 1, B(:,IC), 1)
+       B(1:NROW,IC) = A(1:NROW,IC)
     end do
     RETURN
   end subroutine dacopy

integration/VODE90/cuVODE/source/cuvode_dvjac.F90

@@ -120,8 +120,9 @@ subroutine dvjac(IWM, IERPJ, rwork, vstate)
        CALL JAC (vstate % TN, vstate % Y, 0, 0, &
             rwork % WM(3:3 + VODE_NEQS**2 - 1), VODE_NEQS, vstate % RPAR)
        if (vstate % JSV .EQ. 1) then
-          CALL DCOPYN (LENP, rwork % WM(3:3 + LENP - 1), 1, &
-               rwork % WM(vstate % LOCJS:vstate % LOCJS + LENP - 1), 1)
+          do I = 0, LENP-1
+             rwork % WM(vstate % LOCJS + I) = rwork % WM(3 + I)
+          end do
        endif
     ENDIF
 
@@ -150,16 +151,18 @@ subroutine dvjac(IWM, IERPJ, rwork, vstate)
        vstate % NFE = vstate % NFE + VODE_NEQS
        LENP = VODE_NEQS * VODE_NEQS
        if (vstate % JSV .EQ. 1) then
-          CALL DCOPYN (LENP, rwork % WM(3:3 + LENP - 1), 1, &
-               rwork % WM(vstate % LOCJS:vstate % LOCJS + LENP - 1), 1)
+          do I = 0, LENP-1
+             rwork % WM(vstate % LOCJS + I) = rwork % WM(3 + I)
+          end do
        end if
     ENDIF
 
     IF (JOK .EQ. 1 .AND. (vstate % MITER .EQ. 1 .OR. vstate % MITER .EQ. 2)) THEN
        vstate % JCUR = 0
        LENP = VODE_NEQS * VODE_NEQS
-       CALL DCOPYN (LENP, rwork % WM(vstate % LOCJS:vstate % LOCJS + LENP - 1), 1, &
-            rwork % WM(3:3 + LENP - 1), 1)
+       do I = 0, LENP - 1
+           rwork % WM(3 + I) = rwork % WM(vstate % LOCJS + I)
+       end do
     ENDIF
 
     IF (vstate % MITER .EQ. 1 .OR. vstate % MITER .EQ. 2) THEN

integration/VODE90/cuVODE/source/cuvode_dvnlsd.F90

@@ -137,7 +137,7 @@ subroutine dvnlsd(IWM, NFLAG, rwork, vstate)
     M = 0
     DELP = ZERO
 
-    CALL DCOPYN(VODE_NEQS, rwork % yh(:,1), 1, vstate % Y, 1)
+    vstate % Y(1:VODE_NEQS) = rwork % yh(1:VODE_NEQS,1)
     CALL f_rhs (vstate % TN, vstate % Y, rwork % savf, vstate % RPAR)
     vstate % NFE = vstate % NFE + 1
     IF (vstate % IPUP .LE. 0) GO TO 250
@@ -174,7 +174,7 @@ subroutine dvnlsd(IWM, NFLAG, rwork, vstate)
     do I = 1,VODE_NEQS
        vstate % Y(I) = rwork % YH(I,1) + rwork % SAVF(I)
     end do
-    CALL DCOPYN(VODE_NEQS, rwork % SAVF, 1, rwork % ACOR, 1)
+    rwork % ACOR(1:VODE_NEQS) = rwork % SAVF(1:VODE_NEQS)
 
     GO TO 400
     ! -----------------------------------------------------------------------

integration/VODE90/cuVODE/source/cuvode_dvstep.F90

@@ -293,7 +293,7 @@ subroutine dvstep(IWM, rwork, vstate)
     end do
     vstate % NQWAIT = vstate % NQWAIT - 1
     IF ((vstate % L .EQ. VODE_LMAX) .OR. (vstate % NQWAIT .NE. 1)) GO TO 490
-    CALL DCOPYN(VODE_NEQS, rwork % acor, 1, rwork % yh(:,VODE_LMAX), 1)
+    rwork % yh(1:VODE_NEQS,VODE_LMAX) = rwork % acor(1:VODE_NEQS)
 
     vstate % CONP = vstate % TQ(5)
 490 IF (vstate % ETAMAX .NE. ONE) GO TO 560
@@ -401,7 +401,7 @@ subroutine dvstep(IWM, rwork, vstate)
 620 continue
     vstate % ETA = ETAQP1
     vstate % NEWQ = vstate % NQ + 1
-    CALL DCOPYN(VODE_NEQS, rwork % acor, 1, rwork % yh(:,VODE_LMAX), 1)
+    rwork % yh(1:VODE_NEQS,VODE_LMAX) = rwork % acor(1:VODE_NEQS)
     ! Test tentative new H against THRESH, ETAMAX, and HMXI, then exit. ----
 630 IF (vstate % ETA .LT. THRESH .OR. vstate % ETAMAX .EQ. ONE) GO TO 640
     vstate % ETA = MIN(vstate % ETA,vstate % ETAMAX)