Completed runcvsRoberts_ASAi_dns(). In CVodeB added
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 6 Apr 2018 18:03:48 +0000 (18:03 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 6 Apr 2018 18:03:48 +0000 (18:03 +0000)
CVodeSetMaxNumSteps(tmp_cvB_mem.cv_mem, -1) for unlimited integration steps before call to CVode.  Changed TOUT = 4.0E7 and TB1 = 4.0E7 to TOUT = 4.0E4 and TB1 = 4.0E4 so problem runs in a reasonable amount of time.

git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15445 ba61647d-9d00-f842-95cd-605cb4296b96

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

index 0aba980..ac62889 100644 (file)
@@ -118,19 +118,14 @@ public abstract class CVODES_ASA extends CVODES {
        \r
                /** Problem Constants */\r
                final int NEQ = 3; // Number of equations\r
-               //final double RTOL = 1.0E-6; // 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 RTOL = 1.0E-6; // scalar relative tolerance\r
                final double ATOLl = 1.0E-8; // absolute tolerance for adjoint variables\r
-               final double ATOLq = 1.0E-6; // absolute tolerannce for quuadratures\r
+               final double ATOLq = 1.0E-6; // absolute tolerance for quadratures\r
                final double T0 = 0.0; // initial time\r
-               final double TOUT = 4.0E7; // final time\r
-               final double TB1 = 4.0E7; // starting point for adjoint problem\r
+               //final double TOUT = 4.0E7; // final time\r
+               //final double TB1 = 4.0E7; // starting point for adjoint problem\r
+               final double TOUT = 4.0E4; // final time\r
+           final double TB1 = 4.0E4; // starting point for adjoint problem\r
                final double TB2 = 50.0; // starting point for adjoint problem\r
                final double TBout1 = 40.0; // intermediate t for adjoint problem\r
                final int STEPS = 150; // number of steps between check points\r
@@ -152,7 +147,7 @@ public abstract class CVODES_ASA extends CVODES {
                int ncheck[] = new int[1];\r
                long nst;\r
                int i;\r
-               CVadjCheckPointRec ckpnt[];\r
+               CVadjCheckPointRec ckpnt[] = null;\r
                double reltolB;\r
                double abstolB;\r
                double abstolQB;\r
@@ -360,6 +355,12 @@ public abstract class CVODES_ASA extends CVODES {
                        return;\r
                }\r
                \r
+               // Allow unlimited steps in reaching tout\r
+               flag = CVodeSetMaxNumSteps(cvode_mem, -1);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
                // Call CVodeInitB to allocate internal memory and initialize the\r
                // backward problem.\r
                flag = CVodeInitB(cvode_mem, indexB[0], fB, TB1, yB);\r
@@ -505,6 +506,84 @@ public abstract class CVODES_ASA extends CVODES {
                if (flag != CV_SUCCESS) {\r
                        return;\r
                }\r
+               \r
+               System.out.printf("Backward integration from tB0 = " + TB2 + "\n\n");\r
+               \r
+               // First get results at t = TBout1\r
+               \r
+               flag = CVodeB(cvode_mem, TBout1, CV_NORMAL);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               flag = CVodeGetB(cvode_mem, indexB[0], time, yB);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               flag = CVodeGetAdjY(cvode_mem, TBout1, y);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               System.out.printf("returned t: " + time[0] + "\n");\r
+               System.out.printf("tout: " + TBout1 + "\n");\r
+               System.out.printf("lambda(t): " + yB.data[0] + "  " + yB.data[1] + "  " + yB.data[2] + "\n");\r
+               System.out.printf("y(t): " + y.data[0] + "  " + y.data[1] + "  " + y.data[2] + "\n");\r
+               \r
+               // Then at t = T0\r
+               \r
+               flag = CVodeB(cvode_mem, T0, CV_NORMAL);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               CVodeGetNumSteps(CVodeGetAdjCVodeBmem(cvode_mem, indexB[0]), nstB);\r
+               System.out.printf("Done (nst = " + nstB[0] + ")\n");\r
+               \r
+               flag = CVodeGetB(cvode_mem, indexB[0], time, yB);\r
+               if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+               \r
+               // Call CVodeGetQuadB to get the quadrature solution vector after a\r
+               // successful return from CVodeB.\r
+               flag = CVodeGetQuadB(cvode_mem, indexB[0], time, qB);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               flag = CVodeGetAdjY(cvode_mem, T0, y);\r
+               if (flag < 0) {\r
+                       return;\r
+               }\r
+               \r
+               System.out.printf("returned t: " + time[0] + "\n");\r
+               System.out.printf("lambda(t0): " + yB.data[0] + "  " + yB.data[1] + "  " + yB.data[2] + "\n");\r
+               System.out.printf("y(t0): " + y.data[0] + "  " + y.data[1] + "  " + y.data[2] + "\n");\r
+               System.out.printf("dG/dp: " + (-qB.data[0]) + "  " + (-qB.data[1]) + "  " + (-qB.data[2]) + "\n");\r
+               \r
+               // Free memory\r
+               System.out.printf("Free memory\n\n");\r
+               \r
+               CVodeFree(cvode_mem);\r
+               N_VDestroy(y);\r
+               N_VDestroy(q);\r
+               N_VDestroy(yB);\r
+               N_VDestroy(qB);\r
+               SUNLinSolFree_Dense(LS);\r
+               for (i = 0; i < A.length; i++) {\r
+                       A[i] = null;\r
+               }\r
+               A = null;\r
+               SUNLinSolFree_Dense(LSB);\r
+               for (i = 0; i < AB.length; i++) {\r
+                       AB[i] = null;\r
+               }\r
+               AB = null;\r
+               \r
+               if (ckpnt != null) ckpnt = null;\r
+               data = null;\r
        } \r
        \r
        /*\r
@@ -768,7 +847,7 @@ public abstract class CVODES_ASA extends CVODES {
          case CV_HERMITE:\r
            \r
            ca_mem.ca_IMmalloc = CVAhermiteMalloc_select;\r
-           //ca_mem.ca_IMfree   = CVAhermiteFree;\r
+           ca_mem.ca_IMfree   = CVAhermiteFree_select;\r
            ca_mem.ca_IMget    = CVAhermiteGetY_select;\r
            ca_mem.ca_IMstore  = CVAhermiteStorePnt_select;\r
 \r
@@ -777,7 +856,7 @@ public abstract class CVODES_ASA extends CVODES {
          case CV_POLYNOMIAL:\r
          \r
            ca_mem.ca_IMmalloc = CVApolynomialMalloc_select;\r
-           //ca_mem.ca_IMfree   = CVApolynomialFree;\r
+           ca_mem.ca_IMfree   = CVApolynomialFree_select;\r
            ca_mem.ca_IMget    = CVApolynomialGetY_select;\r
            ca_mem.ca_IMstore  = CVApolynomialStorePnt_select;\r
 \r
@@ -1863,8 +1942,7 @@ public abstract class CVODES_ASA extends CVODES {
          new_cvB_mem.cv_user_data  = null;\r
 \r
          new_cvB_mem.cv_lmem    = null;\r
-         //new_cvB_mem.cv_lfree   = null;\r
-         new_cvB_mem.cv_lfree_select = -1;\r
+         new_cvB_mem.cv_lfree   = -1;\r
          //new_cvB_mem.cv_pmem    = null;\r
          //new_cvB_mem.cv_pfree   = null;\r
 \r
@@ -2082,11 +2160,11 @@ public abstract class CVODES_ASA extends CVODES {
          cvdlsB_mem = new CVDlsMemRecB();\r
 \r
          /* free any existing system solver attached to cvB */\r
-         //if (cvB_mem.cv_lfree_select > 0)  cvB_mem.cv_lfree(cvB_mem);\r
+         //if (cvB_mem.cv_lfree > 0)  cvB_mem.cv_lfree(cvB_mem);\r
 \r
          /* Attach lmemB data and lfreeB function. */\r
          cvB_mem.cv_lmem  = cvdlsB_mem;\r
-         cvB_mem.cv_lfree_select = cvDlsFreeB_select;\r
+         cvB_mem.cv_lfree = cvDlsFreeB_select;\r
 \r
          /* initialize jacB and jacBS pointers */\r
          cvdlsB_mem.jacB  = -1;\r
@@ -2508,6 +2586,7 @@ public abstract class CVODES_ASA extends CVODES {
 \r
                /* Integrate current backward problem */\r
                CVodeSetStopTime(tmp_cvB_mem.cv_mem, ck_mem.ck_t0);\r
+               CVodeSetMaxNumSteps(tmp_cvB_mem.cv_mem, -1);\r
                flag = CVode(tmp_cvB_mem.cv_mem, tBout, tmp_cvB_mem.cv_y, tBret, itaskB);\r
 \r
                /* Set the time at which we will report solution and/or quadratures */\r