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

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

index ca9d371..c711d3c 100644 (file)
@@ -610,7 +610,8 @@ public abstract class CVODES {
     final int cvsRoberts_dnsL = 4;\r
     final int cvsRoberts_FSA_dns = 5;\r
     final int cvsRoberts_ASAi_dns = 6;\r
-    int problem = cvsRoberts_ASAi_dns;\r
+    final int cvsAdvDiff_FSA_non = 7;\r
+    int problem = cvsAdvDiff_FSA_non;\r
     boolean testMode = true;\r
        \r
     // Linear Solver options for runcvsDirectDemo\r
@@ -635,8 +636,8 @@ public abstract class CVODES {
        // Use the following code to call CVODES or CVODES_ASA from another module:\r
        /*boolean testme = true;\r
        Use only one of the following 2 lines\r
-    class CVODEStest extends CVODES { // if not running runcvsRoberts_ASAi_dns()\r
-    class CVODEStest extends CVODES_ASA { // if running runcvsRoberts_ASAi_dns()\r
+    class CVODEStest extends CVODES { // if not running runcvsRoberts_ASAi_dns() or runcvsAdvDiff_FSA_non()\r
+    class CVODEStest extends CVODES_ASA { // if running runcvsRoberts_ASAi_dns() or runcvsAdvDiff_FSA_non()\r
          public CVODEStest() {\r
                  super();\r
          }\r
@@ -2260,6 +2261,36 @@ public abstract class CVODES {
                        ydot[0] = y[1];\r
                        ydot[1] = (ONE - y[0]*y[0]) * P1_ETA * y[1] - y[0];\r
                }\r
+               else if (problem == cvsAdvDiff_FSA_non) {\r
+                       int i;\r
+                       int NEQ = 10;\r
+                   double dx;\r
+                   double ui, ult, urt, hordc, horac, hdiff, hadv;\r
+                   dx = user_data.dx;\r
+                   hordc = user_data.array[0]/(dx*dx);\r
+                   horac = user_data.array[1]/(2.0*dx);\r
+                   \r
+                   /* Loop over all grid points. */\r
+                   for (i=0; i<NEQ; i++) {\r
+\r
+                     /* Extract u at x_i and two neighboring points */\r
+                     ui = y[i];\r
+                     if(i!=0) \r
+                       ult = y[i-1];\r
+                     else\r
+                       ult = ZERO;\r
+                     if(i!=NEQ-1)\r
+                       urt = y[i+1];\r
+                     else\r
+                       urt = ZERO;\r
+\r
+                     /* Set diffusion and advection terms and load into udot */\r
+                     hdiff = hordc*(ult - 2.0*ui + urt);\r
+                     hadv = horac*(urt - ult);\r
+                     ydot[i] = hdiff + hadv;\r
+                   }\r
+\r
+               }\r
                return 0;\r
        }\r
        \r
@@ -2482,6 +2513,7 @@ public abstract class CVODES {
        public class UserData {\r
                CVodeMemRec memRec;\r
                double[] array;\r
+               double dx;\r
        }\r
        \r
          \r
@@ -5062,7 +5094,7 @@ public abstract class CVODES {
     return;\r
   }\r
   \r
-  private double N_VMaxNorm_Serial(NVector x)\r
+  protected double N_VMaxNorm_Serial(NVector x)\r
   {\r
     int i, N;\r
     double max, xd[];\r
@@ -5101,7 +5133,7 @@ public abstract class CVODES {
     return(Math.sqrt(sum/N));\r
   }\r
   \r
-  private void N_VConst_Serial(double c, NVector z)\r
+  protected void N_VConst_Serial(double c, NVector z)\r
   {\r
     int i, N;\r
     double zd[];\r
@@ -11672,7 +11704,7 @@ else                return(snrm);
         * a negative value otherwise.\r
         */\r
     // CVSensRhs1Fn fS1\r
-       private int CVodeSensInit1(CVodeMemRec cv_mem, int Ns, int ism, int fS1_select, NVector yS0[])\r
+       protected int CVodeSensInit1(CVodeMemRec cv_mem, int Ns, int ism, int fS1_select, NVector yS0[])\r
        {\r
          boolean allocOK;\r
          int is;\r
@@ -11925,7 +11957,7 @@ else                return(snrm);
          return(true);\r
        }\r
 \r
-       private int CVodeSensEEtolerances(CVodeMemRec cv_mem)\r
+       protected int CVodeSensEEtolerances(CVodeMemRec cv_mem)\r
        {\r
 \r
          if (cv_mem==null) {\r
@@ -11947,7 +11979,7 @@ else                return(snrm);
          return(CV_SUCCESS);\r
        }\r
 \r
-       private int CVodeSetSensParams(CVodeMemRec cv_mem, double p[], double pbar[], int plist[])\r
+       protected int CVodeSetSensParams(CVodeMemRec cv_mem, double p[], double pbar[], int plist[])\r
        {\r
          int is, Ns;\r
 \r
@@ -12008,7 +12040,7 @@ else                return(snrm);
         * This is just a wrapper that calls CVodeSensDky with k=0.\r
         */\r
         \r
-       private int CVodeGetSens(CVodeMemRec cv_mem, double tret[], NVector ySout[])\r
+       protected int CVodeGetSens(CVodeMemRec cv_mem, double tret[], NVector ySout[])\r
        {\r
          int flag;\r
 \r