Example cvsRoberts_dns_uw works. Added capability for user supplied ewt function...
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 14 Mar 2018 18:59:59 +0000 (18:59 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Wed, 14 Mar 2018 18:59:59 +0000 (18:59 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15415 ba61647d-9d00-f842-95cd-605cb4296b96

mipav/src/gov/nih/mipav/model/algorithms/CVODES.java

index 337f934..8f40ef0 100644 (file)
@@ -552,6 +552,7 @@ public abstract class CVODES {
     \r
     // For cv_mem.cv_efun_select\r
     final int cvEwtSet_select = 1;\r
+    final int cvEwtUser_select1 = 2;\r
     // For cv_mem.cv_linit_select\r
     final int cvDlsInitialize_select = 1;\r
     // For cv_mem.cv_lsolve_select\r
@@ -563,7 +564,8 @@ public abstract class CVODES {
 \r
     final int cvsRoberts_dns = 1;\r
     final int cvSDirectDemo_ls_Problem_1 = 2;\r
-    int problem = cvSDirectDemo_ls_Problem_1;\r
+    final int cvsRoberts_dns_uw = 3;\r
+    int problem = cvsRoberts_dns_uw;\r
     boolean testMode = true;\r
        \r
     // Linear Solver options for runcvsDirectDemo\r
@@ -604,6 +606,10 @@ public abstract class CVODES {
                          NVector tmp2, NVector tmp3) {\r
                  return 0;\r
          }\r
+         \r
+         public int ewt(NVector y, NVector w, CVodeMemRec user_data) {\r
+              return 0;\r
+          }\r
        \r
       }\r
     if (testme) {\r
@@ -636,6 +642,9 @@ public abstract class CVODES {
            if (problem == cvsRoberts_dns) {\r
                runcvsRoberts_dns();\r
            }\r
+           else if (problem == cvsRoberts_dns_uw) {\r
+               runcvsRoberts_dns_uw();\r
+           }\r
            else if (problem == cvSDirectDemo_ls_Problem_1) {\r
                runcvsDirectDemo_Problem_1();\r
            }\r
@@ -912,6 +921,241 @@ public abstract class CVODES {
                return;\r
        }\r
        \r
+       private void runcvsRoberts_dns_uw() {\r
+               /* -----------------------------------------------------------------\r
+                * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh and\r
+                *                Radu Serban @ LLNL\r
+                * -----------------------------------------------------------------\r
+                * Example problem:\r
+                * \r
+                * The following is a simple example problem, with the coding\r
+                * needed for its solution by CVODE. The problem is from\r
+                * chemical kinetics, and consists of the following three rate\r
+                * equations:         \r
+                *    dy1/dt = -.04*y1 + 1.e4*y2*y3\r
+                *    dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*(y2)^2\r
+                *    dy3/dt = 3.e7*(y2)^2\r
+                * on the interval from t = 0.0 to t = 4.e10, with initial\r
+                * conditions: y1 = 1.0, y2 = y3 = 0. The problem is stiff.\r
+                * While integrating the system, we also use the rootfinding\r
+                * feature to find the points at which y1 = 1e-4 or at which\r
+                * y3 = 0.01. This program solves the problem with the BDF method,\r
+                * Newton iteration with the SUNDENSE dense linear solver, and a\r
+                * user-supplied Jacobian routine.\r
+                * It uses a user-supplied function to compute the error weights\r
+                * required for the WRMS norm calculations.\r
+                * Output is printed in decades from t = .4 to t = 4.e10.\r
+                * Run statistics (optional outputs) are printed at the end.\r
+                * -----------------------------------------------------------------*/\r
+               /** Problem Constants */\r
+               final int NEQ = 3; // Number of equations\r
+               final double Y1 = 1.0; // Initial y components\r
+               final double Y2 = 0.0;\r
+               final double Y3 = 0.0;\r
+               //final double RTOL = 1.0E-4; // scalar relative tolerance\r
+               final double RTOL = 1.0E-12;\r
+               //final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
+               final double ATOL1 = 1.0E-12;\r
+               //final double ATOL2 = 1.0E-14;\r
+               final double ATOL2 = 1.0E-15;\r
+               //final double ATOL3 = 1.0E-6;\r
+               final double ATOL3 = 1.0E-12;\r
+               final double T0 = 0.0; // initial time\r
+               final double T1 = 0.4; // first output time\r
+               final double TMULT = 10.0; // output time factor\r
+               //final int NOUT = 12; // number of output times\r
+               final int NOUT = 12;\r
+               int f = cvsRoberts_dns;\r
+               int g = cvsRoberts_dns;\r
+               int Jac = cvsRoberts_dns;\r
+               int ewt_select = cvEwtUser_select1;\r
+               \r
+               // Set initial conditions\r
+               NVector y = new NVector();\r
+               NVector abstol = new NVector();\r
+               double yr[] = new double[]{Y1,Y2,Y3};\r
+               N_VNew_Serial(y, NEQ);\r
+               y.setData(yr);\r
+               double reltol = RTOL; // Set the scalar relative tolerance\r
+               // Set the vector absolute tolerance\r
+               double abstolr[] = new double[]{ATOL1,ATOL2,ATOL3};\r
+               N_VNew_Serial(abstol,NEQ);\r
+               abstol.setData(abstolr);\r
+               CVodeMemRec cvode_mem;\r
+               int flag;\r
+               int flagr;\r
+               double A[][];\r
+               SUNLinearSolver LS;\r
+               double tout;\r
+               int iout;\r
+               double t[] = new double[1];\r
+               int rootsfound[] = new int[2];\r
+               int i;\r
+               \r
+               // Call CVodeCreate to create the solver memory and specify the\r
+               // Backward Differentiation Formula and the use of a Newton\r
+               // iteration\r
+               cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);\r
+               if (cvode_mem == null) {\r
+                   return;     \r
+               }\r
+               // Allow unlimited steps in reaching tout\r
+               flag = CVodeSetMaxNumSteps(cvode_mem, -1);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               // Allow many error test failures\r
+               flag = CVodeSetMaxErrTestFails(cvode_mem, Integer.MAX_VALUE);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeInit to initialize the integrator memory and specify the\r
+               // user's right hand side function in y' = f(t,y), the initial time T0, and\r
+           // the initial dependent variable vector y\r
+               flag = CVodeInit(cvode_mem, f, T0, y);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Use private function to compute error weights\r
+               // CVodeWFtolerances specifies a user-provides function (of type CVEwtFn)\r
+               // which will be called to set the error weight vector.\r
+               flag = CVodeWFtolerances(cvode_mem, ewt_select);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVRootInit to specify the root function g with 2 components\r
+               flag = CVodeRootInit(cvode_mem, 2, g);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Create dense SUNMATRIX for use in linear solver\r
+               // indexed by columns stacked on top of each other\r
+               try {\r
+                   A = new double[NEQ][NEQ];\r
+               }\r
+               catch (OutOfMemoryError e) {\r
+                   MipavUtil.displayError("Out of memory error trying to create A");\r
+                   return;\r
+               }\r
+               \r
+               // Create dense SUNLinearSolver object for use by CVode\r
+               LS = SUNDenseLinearSolver(y, A);\r
+               if (LS == null) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVDlsSetLinearSolver to attach the matrix and linear solver to CVode\r
+               flag = CVDlsSetLinearSolver(cvode_mem, LS, A);\r
+               if (flag != CVDLS_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Set the user-supplied Jacobian routine Jac\r
+               flag = CVDlsSetJacFn(cvode_mem, Jac);\r
+               if (flag != CVDLS_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // In loop, call CVode, print results, and test for error.\r
+               // Break out of loop when NOUT preset output times have been reached.\r
+               System.out.println(" \n3-species kinetics problem\n");\r
+               \r
+               iout = 0;\r
+               tout = T1;\r
+               \r
+               while (true) {\r
+                       flag = CVode(cvode_mem, tout, y, t, CV_NORMAL);\r
+                       switch (iout) {\r
+                       case 0:\r
+                               System.out.println("Example manual: at t = 0.4 y[0] = 0.98517 y[1] = 3.3864E-5 y[2] = 1.4794E-2");\r
+                               break;\r
+                       case 1:\r
+                               System.out.println("Example manual: at t = 4.0 y[0] = 0.90552 y[1] = 2.2405E-5 y[2] = 9.4459E-2");\r
+                               break;\r
+                       case 2:\r
+                               System.out.println("Example manual: at t = 40.0 y[0] = 0.71583 y[1] = 9.1856E-6 y[2] = 0.28416");\r
+                               break;\r
+                       case 3:\r
+                               System.out.println("Example manual: at t = 400.0 y[0] = 0.45052 y[1] = 3.2229E-6 y[2] = 0.54947");\r
+                               break;\r
+                       case 4:\r
+                               System.out.println("Example manual: at t = 4.0E3 y[0] = 0.18317 y[1] = 8.9403E-7 y[2] = 0.81683");\r
+                               break;\r
+                       case 5:\r
+                               System.out.println("Example manual: at t = 4.0E4 y[0] = 3.8977E-2 y[1] = 1.6215E-7 y[2] = 0.96102");\r
+                               break;\r
+                       case 6:\r
+                               System.out.println("Example manual: at t = 4.0E5 y[0] = 4.9387E-3 y[1] = 1.9852E-8 y[2] = 0.99506");\r
+                               break;\r
+                       case 7:\r
+                               System.out.println("Example manual: at t = 4.0E6 y[0] = 5.1684E-4 y[1] = 2.0684E-9 y[2] = 0.99948");\r
+                               break;\r
+                       case 8:\r
+                               System.out.println("Example manual: at t = 4.0E7 y[0] = 5.2039E-5 y[1] = 2.0817E-10 y[2] = 0.99995");\r
+                               break;\r
+                       case 9:\r
+                               System.out.println("Example manual at t = 4.0E8 y[0] = 5.2106E-6 y[1] = 2.0842E-11 y[2] = 0.99999");\r
+                               break;\r
+                       case 10:\r
+                               System.out.println("Example manual at t = 4.0E9 y[0] = 5.1881E-7 y[1] = 2.0752E-12 y[2] = 1.0000");\r
+                               break;\r
+                       case 11:\r
+                               System.out.println("Example manual at t = 4.0E10 y[0] = 6.5181E-8 y[1] = 2.6072E-13 y[2] = 1.0000");\r
+                               break;\r
+                       }\r
+                       System.out.println("MIPAV: at t = " + t[0] + " y[0] = " + y.data[0] + " y[1] = " + y.data[1] + " y[2] = " + y.data[2]);\r
+                       \r
+                       if (flag == CV_ROOT_RETURN) {\r
+                           flagr = CVodeGetRootInfo(cvode_mem, rootsfound)     ;\r
+                           if (flagr == CV_MEM_NULL) {\r
+                               return;\r
+                           }\r
+                           System.out.println("Roots found are " + rootsfound[0] + " and " + rootsfound[1]);\r
+                       }\r
+                       \r
+                       if (flag < 0) {\r
+                               System.out.println("CVode failed with flag = " + flag);\r
+                               break;\r
+                       }\r
+                       \r
+                       if (flag == CV_SUCCESS) {\r
+                               iout++;\r
+                               tout *= TMULT;\r
+                       }\r
+                       \r
+                       if (iout == NOUT) {\r
+                               break;\r
+                       }\r
+               } // while (true)\r
+               \r
+               // Print some final statistics\r
+               PrintFinalStats(cvode_mem);\r
+               \r
+               // Check the solution error\r
+               flag = check_ans(y, reltol, abstol);\r
+               \r
+               // Free y and abstol vectors\r
+               N_VDestroy(y);\r
+               N_VDestroy(abstol);\r
+               \r
+               // Free the integrator memory\r
+               CVodeFree(cvode_mem);\r
+               \r
+               // Free the linear solver memory\r
+               SUNLinSolFree_Dense(LS);\r
+               \r
+               // Free the matrix memory\r
+               for (i = 0; i < NEQ; i++) {\r
+                       A[i] = null;\r
+               }\r
+               A = null;\r
+               return;\r
+       }\r
+       \r
        \r
 \r
        private void PrintFinalStats(CVodeMemRec cv_mem) {\r
@@ -1335,7 +1579,7 @@ public abstract class CVODES {
        private int fTestMode(double t, NVector yv, NVector ydotv, CVodeMemRec user_data) {\r
                double y[] = yv.getData();\r
                double ydot[] = ydotv.getData();\r
-               if (problem == cvsRoberts_dns) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\r
                    ydot[0] = -0.04*y[0] + 1.0E4*y[1]*y[2];\r
                    ydot[2] = 3.0E7*y[1]*y[1];\r
                    ydot[1] = -ydot[0] - ydot[2];\r
@@ -1365,7 +1609,7 @@ public abstract class CVODES {
         */\r
        private int gTestMode(double t, NVector yv, double gout[], CVodeMemRec user_data) {\r
                double y[] = yv.getData();\r
-               if (problem == cvsRoberts_dns) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\r
                    gout[0] = y[0] - 0.0001;\r
                    gout[1] = y[2] - 0.01;\r
                }\r
@@ -1387,7 +1631,7 @@ public abstract class CVODES {
         */\r
        private int JacTestMode(double t, NVector yv, NVector fy, double J[][], CVodeMemRec data, NVector tmp1, NVector tmp2, NVector tmp3) {\r
                double y[] = yv.getData();\r
-               if (problem == cvsRoberts_dns) {\r
+               if ((problem == cvsRoberts_dns) || (problem == cvsRoberts_dns_uw)) {\r
                    J[0][0] = -0.04;\r
                    J[0][1] = 1.0E4 * y[2];\r
                    J[0][2] = 1.0E4 * y[1];\r
@@ -1412,6 +1656,39 @@ public abstract class CVODES {
        \r
        public abstract int Jac(double t, NVector yv, NVector fy, double J[][], CVodeMemRec data, NVector tmp1, NVector tmp2, NVector tmp3);\r
        \r
+       public int ewtTestMode(NVector y, NVector w, CVodeMemRec user_data) {\r
+               if (problem == cvsRoberts_dns_uw) {\r
+                       //final double RTOL = 1.0E-4; // scalar relative tolerance\r
+                       final double RTOL = 1.0E-12;\r
+                       //final double ATOL1 = 1.0E-8; // vector absolute tolerance components\r
+                       final double ATOL1 = 1.0E-12;\r
+                       //final double ATOL2 = 1.0E-14;\r
+                       final double ATOL2 = 1.0E-15;\r
+                       //final double ATOL3 = 1.0E-6;\r
+                       final double ATOL3 = 1.0E-12;\r
+                       int i;\r
+                       double yy, ww, rtol;\r
+                       double atol[] = new double[3];\r
+                       \r
+                       rtol = RTOL;\r
+                       atol[0] = ATOL1;\r
+                       atol[1] = ATOL2;\r
+                       atol[2] = ATOL3;\r
+                       \r
+                       for (i = 0; i < 3; i++) {\r
+                               yy = y.data[i];\r
+                               ww = rtol * Math.abs(yy) + atol[i];\r
+                               if (ww <= 0) {\r
+                                       return -1;\r
+                               }\r
+                               w.data[i] = 1.0/ww;\r
+                       }\r
+               } // if (problem == cvsRoberts_dns_uw)\r
+               return 0;\r
+       }\r
+       \r
+       public abstract int ewt(NVector y, NVector w, CVodeMemRec user_data);\r
+       \r
        \r
        \r
        // Types: struct CVodeMemRec, CVodeMem\r
@@ -2336,6 +2613,31 @@ public abstract class CVODES {
 \r
        return(CV_SUCCESS);\r
      }\r
+     \r
+     private int CVodeWFtolerances(CVodeMemRec cv_mem, int efun_select)\r
+       {\r
+\r
+         if (cv_mem== null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES",\r
+                          "CVodeWFtolerances", MSGCV_NO_MEM);\r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         if (cv_mem.cv_MallocDone == false) {\r
+           cvProcessError(cv_mem, CV_NO_MALLOC, "CVODES",\r
+                          "CVodeWFtolerances", MSGCV_NO_MALLOC);\r
+           return(CV_NO_MALLOC);\r
+         }\r
+\r
+         cv_mem.cv_itol = CV_WF;\r
+\r
+         cv_mem.cv_user_efun = true;\r
+         //cv_mem.cv_efun = efun;\r
+         cv_mem.cv_efun_select = efun_select;\r
+         cv_mem.cv_e_data = null; /* will be set to user_data in InitialSetup */\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
 \r
      private double N_VMin_Serial(NVector x) {\r
         int i;\r
@@ -3103,19 +3405,27 @@ public abstract class CVODES {
       /* Reset and check ewt, ewtQ, ewtS */   \r
       if (cv_mem.cv_nst > 0) {\r
 \r
-        ier = cv_efun(cv_mem.cv_zn[0], cv_mem.cv_ewt, cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
-        if(ier != 0) {\r
-          if (cv_mem.cv_itol == CV_WF)\r
-            cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                           MSGCV_EWT_NOW_FAIL, cv_mem.cv_tn);\r
-          else\r
-            cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
-                           MSGCV_EWT_NOW_BAD, cv_mem.cv_tn);\r
-          istate = CV_ILL_INPUT;\r
-          cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
-          N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
-          break;\r
-        }\r
+          if (testMode && cv_mem.cv_itol == CV_WF) {\r
+                 ier = ewtTestMode(cv_mem.cv_zn[0], cv_mem.cv_ewt, cv_mem.cv_e_data);\r
+          }\r
+          else if ((!testMode) && cv_mem.cv_itol == CV_WF) {\r
+                 ier = ewt(cv_mem.cv_zn[0], cv_mem.cv_ewt, cv_mem.cv_e_data);\r
+          }\r
+          else {\r
+             ier = cv_efun(cv_mem.cv_zn[0], cv_mem.cv_ewt, cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
+          }\r
+          if(ier != 0) {\r
+                 if (cv_mem.cv_itol == CV_WF)\r
+                   cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                                  MSGCV_EWT_NOW_FAIL, cv_mem.cv_tn);\r
+                 else\r
+                   cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVode",\r
+                                  MSGCV_EWT_NOW_BAD, cv_mem.cv_tn);\r
+                 istate = CV_ILL_INPUT;\r
+                 cv_mem.cv_tretlast = tret[0] = cv_mem.cv_tn;\r
+                 N_VScale_Serial(ONE, cv_mem.cv_zn[0], yout);\r
+                 break;\r
+          }\r
 \r
         if (cv_mem.cv_quadr && cv_mem.cv_errconQ) {\r
           ier = cvQuadEwtSet(cv_mem, cv_mem.cv_znQ[0], cv_mem.cv_ewtQ);\r
@@ -3338,8 +3648,18 @@ public abstract class CVODES {
     else                      cv_mem.cv_e_data = cv_mem;\r
 \r
     /* Load initial error weights */\r
-    ier = cv_efun(cv_mem.cv_zn[0], cv_mem.cv_ewt,\r
+    if (testMode && cv_mem.cv_itol == CV_WF) {\r
+        ier = ewtTestMode(cv_mem.cv_zn[0], cv_mem.cv_ewt,\r
+                          cv_mem.cv_e_data);\r
+    }\r
+    else if ((!testMode) && cv_mem.cv_itol == CV_WF) {\r
+       ier = ewt(cv_mem.cv_zn[0], cv_mem.cv_ewt,\r
+                cv_mem.cv_e_data);     \r
+    }\r
+    else {\r
+        ier = cv_efun(cv_mem.cv_zn[0], cv_mem.cv_ewt,\r
                           cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
+    }\r
     if (ier != 0) {\r
       if (cv_mem.cv_itol == CV_WF) \r
         cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "cvInitialSetup",\r
@@ -4535,7 +4855,15 @@ public abstract class CVODES {
    temp2 = cv_mem.cv_acor;\r
 \r
    N_VAbs_Serial(cv_mem.cv_zn[0], temp2);\r
-   cv_efun(cv_mem.cv_zn[0], temp1, cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
+   if (testMode && cv_mem.cv_itol == CV_WF) {\r
+          ewtTestMode(cv_mem.cv_zn[0], temp1, cv_mem.cv_e_data);\r
+   }\r
+   else if ((!testMode) && cv_mem.cv_itol == CV_WF) {\r
+          ewt(cv_mem.cv_zn[0], temp1, cv_mem.cv_e_data);   \r
+   }\r
+   else {\r
+       cv_efun(cv_mem.cv_zn[0], temp1, cv_mem.cv_e_data, cv_mem.cv_efun_select);\r
+   }\r
    N_VInv_Serial(temp1, temp1);\r
    N_VLinearSum_Serial(HUB_FACTOR, temp2, ONE, temp1, temp1);\r
 \r
@@ -10027,5 +10355,4 @@ else                return(snrm);
 \r
 \r
 \r
-\r
 }
\ No newline at end of file