Completed runcvsAdvDiff_FSA_non().
authorilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 6 Apr 2018 21:35:55 +0000 (21:35 +0000)
committerilb@NIH.GOV <ilb@NIH.GOV@ba61647d-9d00-f842-95cd-605cb4296b96>
Fri, 6 Apr 2018 21:35:55 +0000 (21:35 +0000)
git-svn-id: https://citdcbmipav.cit.nih.gov/repos-pub/mipav/trunk@15448 ba61647d-9d00-f842-95cd-605cb4296b96

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

index ac62889..d53de58 100644 (file)
@@ -1,19 +1,25 @@
 package gov.nih.mipav.model.algorithms;\r
 \r
 \r
+import gov.nih.mipav.model.algorithms.CVODES.CVodeMemRec;\r
 import gov.nih.mipav.view.MipavUtil;\r
 \r
 public abstract class CVODES_ASA extends CVODES {\r
        \r
        public CVODES_ASA() {\r
                super();\r
-               runcvsRoberts_ASAi_dns();\r
+               if (problem == cvsRoberts_ASAi_dns) {\r
+                   runcvsRoberts_ASAi_dns();\r
+               }\r
+               else if (problem == cvsAdvDiff_FSA_non) {\r
+                       runcvsAdvDiff_FSA_non();\r
+               }\r
        }\r
        \r
        // In CVODES.java have int problem = cvsRoberts_ASAi_dns;\r
        // Use the following code to call CVODES_ASA from another module:\r
                /*boolean testme = true;\r
-           class CVODEStest extends CVODES_ASA { // if running runcvsRoberts_ASAi_dns()\r
+           class CVODEStest extends CVODES_ASA { // if running runcvsRoberts_ASAi_dns() or runcvsAdvDiff_FSA_non()\r
                  public CVODEStest() {\r
                          super();\r
                  }\r
@@ -61,6 +67,245 @@ public abstract class CVODES_ASA extends CVODES {
                return;\r
            } */\r
        \r
+       private void runcvsAdvDiff_FSA_non() {\r
+               /* -----------------------------------------------------------------\r
+                * Programmer(s): Scott D. Cohen, Alan C. Hindmarsh, George D. Byrne,\r
+                *              and Radu Serban @ LLNL\r
+                * -----------------------------------------------------------------\r
+                * Example problem:\r
+                *\r
+                * The following is a simple example problem, with the program for\r
+                * its solution by CVODES. The problem is the semi-discrete form of\r
+                * the advection-diffusion equation in 1-D:\r
+                *   du/dt = q1 * d^2 u / dx^2 + q2 * du/dx\r
+                * on the interval 0 <= x <= 2, and the time interval 0 <= t <= 5.\r
+                * Homogeneous Dirichlet boundary conditions are posed, and the\r
+                * initial condition is:\r
+                *   u(x,y,t=0) = x(2-x)exp(2x).\r
+                * The PDE is discretized on a uniform grid of size MX+2 with\r
+                * central differencing, and with boundary values eliminated,\r
+                * leaving an ODE system of size NEQ = MX.\r
+                * This program solves the problem with the option for nonstiff\r
+                * systems: ADAMS method and functional iteration.\r
+                * It uses scalar relative and absolute tolerances.\r
+                * Output is printed at t = .5, 1.0, ..., 5.\r
+                * Run statistics (optional outputs) are printed at the end.\r
+                *\r
+                * Optionally, CVODES can compute sensitivities with respect to the\r
+                * problem parameters q1 and q2.\r
+                * Any of three sensitivity methods (SIMULTANEOUS, STAGGERED, and\r
+                * STAGGERED1) can be used and sensitivities may be included in the\r
+                * error test or not (error control set on FULL or PARTIAL,\r
+                * respectively).\r
+                *\r
+                * Execution:\r
+                *\r
+                * If no sensitivities are desired:\r
+                *    % cvsAdvDiff_FSA_non -nosensi\r
+                * If sensitivities are to be computed:\r
+                *    % cvsAdvDiff_FSA_non -sensi sensi_meth err_con\r
+                * where sensi_meth is one of {sim, stg, stg1} and err_con is one of\r
+                * {t, f}.\r
+                * -----------------------------------------------------------------*/\r
+               \r
+               /* Problem Constants */\r
+               double XMAX = 2.0;   /* domain boundary           */\r
+               int MX = 10;   /* mesh dimension            */\r
+               int NEQ = MX;  /* number of equations       */\r
+               double ATOL = 1.e-5; /* scalar absolute tolerance */\r
+               double T0 = 0.0;   /* initial time              */\r
+               double T1 = 0.5;   /* first output time         */\r
+               double DTOUT = 0.5;   /* output time increment     */\r
+               int NOUT = 10;    /* number of output times    */\r
+\r
+               int NP = 2;\r
+               int NS = 2;\r
+\r
+               CVodeMemRec cvode_mem;\r
+               UserData data;\r
+               double dx, reltol, abstol, tout;\r
+               double t[] = new double[1];\r
+           NVector u;\r
+               int iout, flag;\r
+\r
+               double pbar[];\r
+               int is, plist[];\r
+               NVector uS[];\r
+               double udata[];\r
+               int i;\r
+               double x;\r
+               int f = cvsAdvDiff_FSA_non;\r
+               int qu;\r
+               double hu;\r
+\r
+               cvode_mem = null;\r
+               data = null;\r
+               u = null;\r
+               pbar = null;\r
+               plist = null;\r
+               uS = null;\r
+               \r
+               // In CVODES.java set sensi to true or false.  If sensi is true, set sensi_meth to\r
+               // CV_SIMULTANEOUS, CV_STAGGERED, or CV_STAGGERED1.  Set err_con to true or false.\r
+               data = new UserData();\r
+               data.array = new double[NP];\r
+        dx = data.dx = XMAX/((double)(MX+1));\r
+        data.array[0] = 1.0;\r
+        data.array[1] = 0.5;\r
+        \r
+        // Allocate and set initial states\r
+        u = new NVector();\r
+        N_VNew_Serial(u,NEQ);\r
+        \r
+        /* Set pointer to data array and get local length of u. */\r
+           udata = u.data;\r
+\r
+           /* Load initial profile into u vector */\r
+           for (i=0; i<NEQ; i++) {\r
+             x = (i+1)*dx;\r
+             udata[i] = x*(XMAX - x)*Math.exp(2.0*x);\r
+           }  \r
+           \r
+           /* Set integration tolerances */\r
+           reltol = ZERO;\r
+           abstol = ATOL;\r
+\r
+           /* Create CVODES object */\r
+           cvode_mem = CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);\r
+           if (cvode_mem == null) {\r
+                   return;     \r
+        }\r
+           \r
+           cvode_mem.cv_user_data = data;\r
+           \r
+           /* Allocate CVODES memory */\r
+           flag = CVodeInit(cvode_mem, f, T0, u);\r
+           if (flag != CV_SUCCESS) {\r
+                       return;\r
+               }\r
+           \r
+           flag = CVodeSStolerances(cvode_mem, reltol, abstol);\r
+           if (flag < 0) {\r
+                 return;\r
+           }\r
+           \r
+           System.out.printf("\n1-D advection-diffusion equation, mesh size = " + MX + "\n");\r
+\r
+           /* Sensitivity-related settings */\r
+           if(sensi) {\r
+\r
+             plist = new int[NS];\r
+             for(is=0; is<NS; is++) plist[is] = is;\r
+\r
+             pbar  = new double[NS];\r
+             for(is=0; is<NS; is++) pbar[is] = data.array[plist[is]];\r
+\r
+             uS = N_VCloneVectorArray_Serial(NS, u);\r
+             if (uS == null) {\r
+                 return;\r
+             }\r
+             for(is=0;is<NS;is++)\r
+                 N_VConst_Serial(ZERO, uS[is]);\r
+\r
+             flag = CVodeSensInit1(cvode_mem, NS, sensi_meth, -1, uS);\r
+             if (flag < 0) {\r
+                         return;\r
+                 }\r
+\r
+             flag = CVodeSensEEtolerances(cvode_mem);\r
+             if (flag < 0) {\r
+                         return;\r
+                 }\r
+\r
+             // Set sensitivity analysis optional inputs.\r
+             // Call CVodeSetSensErrCon to specify the error control strategy for\r
+             // sensitivity variables.\r
+                 cvode_mem.cv_errconS = err_con;\r
+\r
+             flag = CVodeSetSensDQMethod(cvode_mem, CV_CENTERED, ZERO);\r
+             if (flag != CV_SUCCESS) {\r
+                         return;\r
+                 }\r
+\r
+             flag = CVodeSetSensParams(cvode_mem, data.array, pbar, plist);\r
+             if (flag < 0) {\r
+                         return;\r
+                 }\r
+\r
+             System.out.printf("Sensitivity: YES ");\r
+             if(sensi_meth == CV_SIMULTANEOUS)   \r
+                 System.out.printf("( SIMULTANEOUS +");\r
+             else \r
+               if(sensi_meth == CV_STAGGERED) System.out.printf("( STAGGERED +");\r
+               else                           System.out.printf("( STAGGERED1 +");   \r
+             if(err_con) System.out.printf(" FULL ERROR CONTROL )");\r
+             else        System.out.printf(" PARTIAL ERROR CONTROL )");\r
+\r
+         } else {\r
+\r
+             System.out.printf("Sensitivity: NO ");\r
+\r
+         }\r
+           \r
+         /* In loop over output points, call CVode, print results, test for error */\r
+\r
+         \r
+         for (iout=1, tout=T1; iout <= NOUT; iout++, tout += DTOUT) {\r
+           flag = CVode(cvode_mem, tout, u, t, CV_NORMAL);\r
+           if(flag < 0) break;\r
+           System.out.println("t = " + t[0]);\r
+           // Returns the order on the last successful step\r
+           qu = cvode_mem.cv_qu;\r
+           // Return the step size used on the last successful step\r
+           hu = cvode_mem.cv_hu;\r
+           System.out.println("Step order qu = " + qu);\r
+           System.out.println("Step size hu = " + hu);\r
+           System.out.println("Number of integrations steps = " + cvode_mem.cv_nst);\r
+           System.out.println("Max norm of solution = " + N_VMaxNorm_Serial(u));\r
+           if (sensi) {\r
+             flag = CVodeGetSens(cvode_mem, t, uS);\r
+             if(flag < 0) break;\r
+             System.out.println("Max norm of sensitivity 1 -= " + N_VMaxNorm_Serial(uS[0]));\r
+             System.out.println("Max norm of sensitivity 2 -= " + N_VMaxNorm_Serial(uS[1]));\r
+           } \r
+           System.out.printf("------------------------------------------------------------\n");\r
+         }\r
+\r
+         /* Print final statistics */\r
+         PrintFinalStats(cvode_mem, sensi);\r
+\r
+         /* Free memory */\r
+         N_VDestroy(u);\r
+         if (sensi) {\r
+           N_VDestroyVectorArray_Serial(uS, NS);\r
+           plist = null;\r
+           pbar = null;\r
+         }\r
+         data.array = null;\r
+         data = null;\r
+         CVodeFree(cvode_mem);\r
+       }\r
+       \r
+       private void PrintFinalStats(CVodeMemRec cv_mem, boolean sensi) {\r
+               System.out.println("Final statistics:");\r
+               System.out.println("Number of integrations steps = " + cv_mem.cv_nst);\r
+               System.out.println("Number of calls to f = " + cv_mem.cv_nfe);\r
+               System.out.println("Number of calls to the linear solver setup routine = " + cv_mem.cv_nsetups);\r
+               System.out.println("Number of error test failures = " + cv_mem.cv_netf[0]);\r
+           System.out.println("Number of nonlinear solver iterations = " + cv_mem.cv_nni);\r
+           System.out.println("Number of nonlinear solver convergence failures = " + cv_mem.cv_ncfn[0]);\r
+           if (sensi) {\r
+               System.out.println("Number of calls to the sensitivity right hand side routine = " + cv_mem.cv_nfSe);\r
+               System.out.println("Number of calls to the user f routine due to finite difference evaluations ");\r
+               System.out.println("of the sensitivity equations = " + cv_mem.cv_nfeS);\r
+               System.out.println("Number of calls made to the linear solver's setup routine due to ");\r
+            System.out.println("sensitivity computations = " + cv_mem.cv_nsetupsS);\r
+            System.out.println("Number of local error test failures for sensitivity variables = " + cv_mem.cv_netfS[0]);\r
+            System.out.println("Total number of nonlinear iterations for sensitivity variables = " + cv_mem.cv_nniS);\r
+            System.out.println("Total number of nonlinear convergence failures for sensitivity variables = " + cv_mem.cv_ncfnS[0]);\r
+           }\r
+       }\r
+       \r
        private void runcvsRoberts_ASAi_dns() {\r
                /* -----------------------------------------------------------------\r
                 * Programmer(s): Radu Serban @ LLNL\r
@@ -585,6 +830,7 @@ public abstract class CVODES_ASA extends CVODES {
                if (ckpnt != null) ckpnt = null;\r
                data = null;\r
        } \r
+\r
        \r
        /*\r
         * CVodeQuadInit\r
@@ -3418,5 +3664,30 @@ public abstract class CVODES_ASA extends CVODES {
 \r
          return(CV_SUCCESS);\r
        }\r
+       \r
+       private int CVodeSetSensDQMethod(CVodeMemRec cv_mem, int DQtype, double DQrhomax)\r
+       {\r
+\r
+         if (cv_mem==null) {\r
+           cvProcessError(null, CV_MEM_NULL, "CVODES", "CVodeSetSensDQMethod", MSGCV_NO_MEM);    \r
+           return(CV_MEM_NULL);\r
+         }\r
+\r
+         if ( (DQtype != CV_CENTERED) && (DQtype != CV_FORWARD) ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSetSensDQMethod", MSGCV_BAD_DQTYPE);    \r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         if (DQrhomax < ZERO ) {\r
+           cvProcessError(cv_mem, CV_ILL_INPUT, "CVODES", "CVodeSetSensDQMethod", MSGCV_BAD_DQRHO);    \r
+           return(CV_ILL_INPUT);\r
+         }\r
+\r
+         cv_mem.cv_DQtype = DQtype;\r
+         cv_mem.cv_DQrhomax = DQrhomax;\r
+\r
+         return(CV_SUCCESS);\r
+       }\r
+\r
 \r
 }
\ No newline at end of file