Matthew Turk avatar Matthew Turk committed 9ce3fe2

Moving allocations outside; these will get blown away by new versions of Dengo.

Comments (0)

Files changed (4)

src/enzo/BE_chem_solve.C

 int BE_chem_solve(rhs_f f, jac_f J,
 		  double *u, double dt, double *rtol, 
                   double *atol, int nstrip, int nchem, 
-		  double *scaling, void *sdata) {
+		  double *scaling, void *sdata,
+          double *u0, double *s, double *gu, double *Ju) {
 
   // local variables
   int i, j, ix, isweep, ier, ONE=1, ioff;
   //fprintf(stderr, "nchem = %d, nstrip = %d\n", nchem, nstrip);
 
   // create/initialize temporary arrays
-  double *u0 = new double[nchem*nstrip];        // initial state
-  double *s  = new double[nchem];               // Newton update (each cell)
-  double *gu = new double[nchem*nstrip];        // nonlinear residual
-  double *Ju = new double[nchem*nchem*nstrip];  // Jacobian
+  //double *u0 = new double[nchem*nstrip];        // initial state
+  //double *s  = new double[nchem];               // Newton update (each cell)
+  //double *gu = new double[nchem*nstrip];        // nonlinear residual
+  //double *Ju = new double[nchem*nchem*nstrip];  // Jacobian
 
   //fprintf(stderr, "Made it here.\n");
   
       ioff = ix*nchem;
 
       // solve for Newton update
-      if (Gauss_Elim(&(Ju[ix*nchem*nchem]), s, &(gu[ioff]), nchem) != 0)
-    return 1;
+       if (Gauss_Elim(&(Ju[ix*nchem*nchem]), s, &(gu[ioff]), nchem) != 0) {
+         unsolved = 1;
+         break;
+       }
 	/*ENZO_FAIL("Error in Gauss_Elim");*/
 
       // update solution in this cell
 	  printf("BE_chem_solve ERROR: NaN in iteration %i (cell %i, species %i)\n",isweep,ix,i);
 	  // rescale back to input variables
 	  for (i=0; i<nstrip*nchem; i++)  u[i] *= scaling[i];
-	  return 1;
+	  unsolved = 1;
+      break;
 	}
       } // i loop
     } // ix loop
   } // end newton iterations
 
   // free temporary arrays
-  delete[] u0;
-  delete[] s;
-  delete[] gu;
-  delete[] Ju;
+  //delete[] u0;
+  //delete[] s;
+  //delete[] gu;
+  //delete[] Ju;
 
   // rescale back to input variables
   for (i=0; i<nstrip*nchem; i++)  u[i] *= scaling[i];

src/enzo/Grid_SolveDengo.C

     }
 
     // call dengo on this strip
-    fprintf(stderr, "Calling dengo_evolve now.\n");
+    /*fprintf(stderr, "Calling dengo_evolve now.\n");*/
     ttot = dengo_evolve_oxygen(dtFixed*TimeUnits, DengoData.dt, strip, rtol, atol,
 			       strip_len, (oxygen_data *) DengoData.dengo_solver_data);
     if (ttot < dtFixed*TimeUnits) // Did we make it?

src/enzo/oxygen_solver.C

 */
 
 
+#include <alloca.h>
 #include "oxygen_solver.h"
 
 oxygen_data *oxygen_setup_data(void) {
             double *rtol, double *atol, long long dims, oxygen_data *data) {
     int i, j;
     hid_t file_id;
-    fprintf(stderr, "  ncells = % 3i\n", (int) dims);
+    /*fprintf(stderr, "  ncells = % 3i\n", (int) dims);*/
 
     int N = 11;
     rhs_f f = calculate_rhs_oxygen;
     double *prev = (double *) alloca(dims * N * sizeof(double));
     for (i = 0; i < dims * N; i++) scale[i] = 1.0;
     for (i = 0; i < dims * N; i++) prev[i] = input[i];
+    double *u0 = (double *) alloca(N*dims*sizeof(double));
+    double *s  = (double *) alloca(N*sizeof(double));
+    double *gu = (double *) alloca(N*dims*sizeof(double));
+    double *Ju = (double *) alloca(N*N*dims*sizeof(double));
     while (ttot < dtf) {
-        int rv = BE_chem_solve(f, jf, input, dt, rtol, atol, dims, N, scale, (void *) data);
+        int rv = BE_chem_solve(f, jf, input, dt, rtol, atol, dims, N, scale, (void *) data,
+            u0, s, gu, Ju);
 	/*
 	fprintf(stderr, "Made it passed BE_chem_solve.\n");
       
         }
         niter++;
     }
-    fprintf(stderr, "End: %0.5g / %0.5g (%0.5g)\n",
-        ttot, dtf, dtf-ttot);
+    /*fprintf(stderr, "End: %0.5g / %0.5g (%0.5g)\n",
+        ttot, dtf, dtf-ttot);*/
 
     return ttot;
 }

src/enzo/oxygen_solver.h

 typedef int(*rhs_f)(double *, double *, int, int, void *);
 typedef int(*jac_f)(double *, double *, int, int, void *);
 int BE_chem_solve(rhs_f f, jac_f J, double *u, double dt, double *rtol, 
-                  double *atol, int nstrip, int nchem, double *scaling, void *sdata);
+                  double *atol, int nstrip, int nchem, double *scaling, void *sdata,
+                  double *, double *, double *, double *);
 
 
 
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.