Commits

Devin Silvia  committed c8300f1

Fixing cell iteration and trying to get Collapse Test working with Dengo.

Made changes to the way we break a grid up into strips and added code to
try to get CollapseTest working. Currently throws a malloc error in
BE_chem_solve.C after a few interations through multiple strips.

  • Participants
  • Parent commits c669380

Comments (0)

Files changed (5)

File src/enzo/BE_chem_solve.C

   // rescale input to normalized variables
   for (i=0; i<nstrip*nchem; i++)  u[i] /= scaling[i];
 
+  //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 *Ju = new double[nchem*nchem*nstrip];  // Jacobian
+
+  //fprintf(stderr, "Made it here.\n");
   
   for (i=0; i<nstrip*nchem; i++)        u0[i] = u[i];
   for (i=0; i<nstrip*nchem; i++)        gu[i] = 0.0;
       // check error in this cell (max norm)
       for (i=0; i<nchem; i++) {
 	if ( fabs(s[i]) > (atol[ioff+i] + rtol[ioff+i] * fabs(u[ioff+i]))) {
-      /*fprintf(stderr, "Unsolved[%d]: %d % 0.5g % 0.5g\n",
-         ix, i, s[i], atol[ioff+i] + rtol[ioff+i] * fabs(u[ioff+i]));*/
+	  /*fprintf(stderr, "Unsolved[%d]: %d % 0.5g % 0.5g\n",
+	    ix, i, s[i], atol[ioff+i] + rtol[ioff+i] * fabs(u[ioff+i]));*/
 	  unsolved = 1;
 	  break;
 	}

File src/enzo/CollapseTestInitialize.C

 void AddLevel(LevelHierarchyEntry *Array[], HierarchyEntry *Grid, int level);
 int RebuildHierarchy(TopGridData *MetaData,
 		     LevelHierarchyEntry *LevelArray[], int level);
+int InitializeDengoData(grid *Grid);
 
 int CollapseTestInitialize(FILE *fptr, FILE *Outfptr, 
 			  HierarchyEntry &TopGrid, TopGridData &MetaData)
   const char *HDIName   = "HDI_Density";
   const char *MetalName = "Metal_Density";
 
+  // Dengo species names
+  const char *OIName    = "OI_Density";
+  const char *OIIName    = "OII_Density";
+  const char *OIIIName    = "OIII_Density";
+  const char *OIVName    = "OIV_Density";
+  const char *OVName    = "OV_Density";
+  const char *OVIName    = "OVI_Density";
+  const char *OVIIName    = "OVII_Density";
+  const char *OVIIIName    = "OVIII_Density";
+  const char *OIXName    = "OIX_Density";
+
   /* declarations */
 
   char  line[MAX_LINE_LENGTH];
     ENZO_FAIL("Error in CollapseTestInitializeGrid.");
   }
 
+  /* Initialize Dengo Data */
+  if (DengoChemistryModel)
+    InitializeDengoData(TopGrid.GridData);
+
   /* Convert minimum initial overdensity for refinement to mass
      (unless MinimumMass itself was actually set). */
 
   if (CollapseTestUseMetals)
     DataLabel[count++] = (char*) MetalName;
 
+  if (DengoChemistryModel) {
+    if (MultiSpecies == 0)
+      DataLabel[count++] = (char*) ElectronName;
+    DataLabel[count++] = (char*) OIName;
+    DataLabel[count++] = (char*) OIIName;
+    DataLabel[count++] = (char*) OIIIName;
+    DataLabel[count++] = (char*) OIVName;
+    DataLabel[count++] = (char*) OVName;
+    DataLabel[count++] = (char*) OVIName;
+    DataLabel[count++] = (char*) OVIIName;
+    DataLabel[count++] = (char*) OVIIIName;
+    DataLabel[count++] = (char*) OIXName;
+  } //if DengoChemistryModel
+
   for (i = 0; i < count; i++)
     DataUnits[i] = NULL;
 

File src/enzo/Grid_CollapseTestInitializeGrid.C

   int DeNum, HINum, HIINum, HeINum, HeIINum, HeIIINum, HMNum, H2INum, H2IINum,
     DINum, DIINum, HDINum, MetalNum;
 
+  // Dengo
+  int OINum, OIINum, OIIINum, OIVNum, OVNum,
+    OVINum, OVIINum, OVIIINum, OIXNum;
+  float ion_fraction = 1e-8; //temporary solution
+
   /* create fields */
 
   NumberOfBaryonFields = 0;
   if (SphereUseColour)
     FieldType[NumberOfBaryonFields++] = Metallicity; /* fake it with metals */
 
+  // Dengo Species -- at some point, we might need to figure out
+  // how to template this part
+  if (DengoChemistryModel) {
+    if (MultiSpecies == 0)
+      FieldType[DeNum   = NumberOfBaryonFields++] = ElectronDensity;
+    FieldType[OINum     = NumberOfBaryonFields++] = OIDensity;
+    FieldType[OIINum    = NumberOfBaryonFields++] = OIIDensity;
+    FieldType[OIIINum   = NumberOfBaryonFields++] = OIIIDensity;
+    FieldType[OIVNum    = NumberOfBaryonFields++] = OIVDensity;
+    FieldType[OVNum     = NumberOfBaryonFields++] = OVDensity;
+    FieldType[OVINum    = NumberOfBaryonFields++] = OVIDensity;
+    FieldType[OVIINum   = NumberOfBaryonFields++] = OVIIDensity;
+    FieldType[OVIIINum  = NumberOfBaryonFields++] = OVIIIDensity;
+    FieldType[OIXNum    = NumberOfBaryonFields++] = OIXDensity;
+  }
+
   /* Return if this doesn't concern us. */
 
   if (ProcessorNumber != MyProcessorNumber) {
 	    }
 	  }
 
+	  // Set the Dengo Species values ... going to hard code these number for now
+	  // we might want to do something along the lines of "TestProblemData.HI_Fraction" type stuff
+	  // Need to think about whether or not a factor of 16 needs to be included...
+	  // Also, we're not going to tie this into any sort of mass fraction of oxygen or metallicity
+	  // for now.
+	  if (DengoChemistryModel) {
+	    BaryonField[OINum][n]    = BaryonField[0][n];
+	    BaryonField[OIINum][n]   = ion_fraction * BaryonField[0][n];
+	    BaryonField[OIIINum][n]  = ion_fraction * BaryonField[0][n];
+	    BaryonField[OIVNum][n]   = ion_fraction * BaryonField[0][n];
+	    BaryonField[OVNum][n]    = ion_fraction * BaryonField[0][n];
+	    BaryonField[OVINum][n]   = ion_fraction * BaryonField[0][n];
+	    BaryonField[OVIINum][n]  = ion_fraction * BaryonField[0][n];
+	    BaryonField[OVIIINum][n] = ion_fraction * BaryonField[0][n];
+	    BaryonField[OIXNum][n]   = ion_fraction * BaryonField[0][n];
+	    // Adjust the OI value using species conservation
+	    BaryonField[OINum][n] -= (BaryonField[OIINum][n] + BaryonField[OIIINum][n] + \
+				      BaryonField[OIVNum][n] + BaryonField[OVNum][n] + \
+				      BaryonField[OVINum][n] + BaryonField[OVIINum][n] + \
+				      BaryonField[OVIIINum][n] + BaryonField[OIXNum][n]);
+	    // Adjust the electron density contributions from oxygen
+	    if (MultiSpecies == 0)
+	      BaryonField[DeNum][n] = (BaryonField[OIINum][n] + (2.0 * BaryonField[OIIINum][n]) + \
+				       (3.0 * BaryonField[OIVNum][n]) + (4.0 * BaryonField[OVNum][n]) + \
+				       (5.0 * BaryonField[OVINum][n]) + (6.0 * BaryonField[OVIINum][n]) + \
+				       (7.0 * BaryonField[OVIIINum][n]) + (8.0 * BaryonField[OIXNum][n]));
+	    BaryonField[DeNum][n] += (BaryonField[OIINum][n] + (2.0 * BaryonField[OIIINum][n]) + \
+				      (3.0 * BaryonField[OIVNum][n]) + (4.0 * BaryonField[OVNum][n]) + \
+				      (5.0 * BaryonField[OVINum][n]) + (6.0 * BaryonField[OVIINum][n]) + \
+				      (7.0 * BaryonField[OVIIINum][n]) + (8.0 * BaryonField[OIXNum][n]));
+	    
+	  } //if (DengoChemistryModel)
+
 	  /* If there are metals, set it. */
 
 	  if (SphereUseMetals)

File src/enzo/Grid_SolveDengo.C

   for (dim = 0; dim < GridRank; dim++)
     size *= GridDimension[dim];
 
+  fprintf(stderr, "\n size = %"ISYM"\n", size);
+
   // get units 
   FLOAT a=1.0, dadt;
   float TemperatureUnits=1.0, DensityUnits=1.0, LengthUnits=1.0, 
 
   //--------------
   // call dengo in strips of spatial length 1024
-  int icur=0, ierr=0, strip_len;
+  int icur, ierr=0, strip_len;
   float *strip = (float *) malloc(1024 * DengoData.NumberOfFields * sizeof(float));
   float *atol = (float *) malloc(1024 * DengoData.NumberOfFields * sizeof(float));
   float *rtol = (float *) malloc(1024 * DengoData.NumberOfFields * sizeof(float));
-  DengoData.dt = -1.0;
+  double ttot;
   for (icur=0; icur<size; icur+=1024) {
+    DengoData.dt = -1.0;
+    ttot = 0.0;
 
     // determine strip length
     strip_len = (1024 < size-icur) ? 1024 : size-icur;
-    printf("strip_len = %"ISYM"\n", strip_len);
+    //fprintf(stderr, "strip_len = %"ISYM"\n", strip_len);
 
     // fill strip, atol, and rtol in order specified by Dengo's "Field" list
     for (ifield=0; ifield<DengoData.NumberOfFields; ifield++) {
-      for (i=0; i<strip_len; i++) {
+      for (i=icur; i < icur+strip_len; i++) { 
 	strip[i*DengoData.NumberOfFields + ifield] = BaryonField[DengoData.Field[ifield]][i];
 	atol[i*DengoData.NumberOfFields + ifield] = BaryonField[DengoData.Field[ifield]][i]
 	                                            * DengoData.ToleranceFactor;
 	rtol[i*DengoData.NumberOfFields + ifield] = DengoData.ToleranceFactor;
+	/*
+	fprintf(stderr, "i = %d, icur = %d, icur+strip_len = %d, strip value = %0.5g\n", i, icur, icur+strip_len,
+		strip[i*DengoData.NumberOfFields + ifield]);
+	*/
       }
     }
 
     // call dengo on this strip
-    double ttot;
     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);
 
     // unpack Dengo solution back into BaryonField arrays
     for (ifield=0; ifield<DengoData.NumberOfFields; ifield++) {
-      for (i=0; i<strip_len; i++) 
+      for (i=icur; i < icur+strip_len; i++) 
 	BaryonField[DengoData.Field[ifield]][i] = strip[i*DengoData.NumberOfFields + ifield];
     }
 

File src/enzo/oxygen_solver.C

     for (i = 0; i < dims * N; i++) prev[i] = input[i];
     while (ttot < dtf) {
         int rv = BE_chem_solve(f, jf, input, dt, rtol, atol, dims, N, scale, (void *) data);
-        /*
+	/*
+	fprintf(stderr, "Made it passed BE_chem_solve.\n");
+      
         fprintf(stderr, "Return value [%d]: %i.  %0.5g / %0.5g = %0.5g (%0.5g)\n",
                 niter, rv, ttot, dtf, ttot/dtf, dt);
         fprintf(stderr, "Value[80] = %0.5g %0.5g %0.5g\n",
         for (i = 0; i < dims * N; i++) {
             if (input[i] < 0) {
                 rv = 1;
+		//fprintf(stderr, "rv = 1, input[%d] = %0.5g\n", i, input[i]);
                 break;
             }
         }