Commits

Daniel Reynolds committed ab2f993

Added Grid_SolveDengo routine; still needs work to flesh out API for Dengo data structure and solver interface.

Comments (0)

Files changed (2)

 
    int SolveRateAndCoolEquations(int RTCoupledSolverIntermediateStep);
 
+/* Solve the joint rate and radiative cooling/heating equations via Dengo */
+
+   int SolveDengo();
+
 /* Compute densities of various species for RadiationFieldUpdate. */
 
    int RadiationComputeDensities(int level);

src/enzo/Grid_SolveDengo.C

+/***********************************************************************
+/
+/  GRID CLASS (SOLVE THE COOLING/HEATING AND RATE EQUATIONS VIA DENGO)
+/
+/  written by: Dan Reynolds, Devin Silva & Matt Turk
+/  date:       December 2012
+/  modified1:  
+/
+/  PURPOSE:
+/
+/  RETURNS:
+/    SUCCESS or FAIL
+/
+************************************************************************/
+
+#include <stdio.h>
+#include <math.h>
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "Grid.h"
+#include "fortran.def"
+#include "CosmologyParameters.h"
+#include "Gadget.h"
+
+// function prototypes
+int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
+int GetUnits(float *DensityUnits, float *LengthUnits,
+	     float *TemperatureUnits, float *TimeUnits,
+	     float *VelocityUnits, FLOAT Time);
+int RadiationFieldCalculateRates(FLOAT Time);
+int FindField(int field, int farray[], int numfields);
+
+
+int grid::SolveDengo()
+{
+  // return if this doesn't concern us.
+  if (!(MultiSpecies && RadiativeCooling)) return SUCCESS;
+  if (ProcessorNumber != MyProcessorNumber)
+    return SUCCESS;
+  if (NumberOfBaryonFields == 0)
+    return SUCCESS;
+
+  this->DebugCheck("SolveDengo");
+
+  // declarations
+  int IENum;
+  FLOAT a = 1.0, dadt;
+    
+  // find internal energy density field
+  if ((GENum = FindField(InternalEnergy, FieldType, NumberOfBaryonFields)) < 0)
+    ENZO_FAIL("Could not find internal energy field.\n");
+
+  // find photo-ionization fields
+  int kphHINum, kphHeINum, kphHeIINum, kdissH2INum;
+  int gammaNum;
+  IdentifyRadiativeTransferFields(kphHINum, gammaNum, kphHeINum, 
+				  kphHeIINum, kdissH2INum);
+
+  // calculate the rates due to the radiation field
+  if (!GadgetEquilibriumCooling) {
+    if (RadiationFieldCalculateRates(Time+0.5*dtFixed) == FAIL)
+      ENZO_FAIL("Error in RadiationFieldCalculateRates.");
+  }
+
+  // compute size of the current grid
+  int i, dim, size = 1;
+  for (dim = 0; dim < GridRank; dim++)
+    size *= GridDimension[dim];
+
+  // get units 
+  float TemperatureUnits = 1.0, DensityUnits = 1.0, LengthUnits = 1.0, 
+    VelocityUnits = 1.0, TimeUnits = 1.0, aUnits = 1.0;
+  if (ComovingCoordinates) {
+    if (CosmologyComputeExpansionFactor(Time+0.5*dtFixed, &a, &dadt) == FAIL)
+      ENZO_FAIL("Error in CosmologyComputeExpansionFactors.");
+    aUnits = 1.0/(1.0 + InitialRedshift);
+  }
+  if (GetUnits(&DensityUnits, &LengthUnits, &TemperatureUnits,
+	       &TimeUnits, &VelocityUnits, Time) == FAIL) {
+    ENZO_FAIL("Error in GetUnits.");
+  }
+
+  //--------------
+  // rescale all relevant fields to CGS
+  //    internal energy
+  float rscale = LengthUnits*LengthUnits/TimeUnits/TimeUnits;
+  for (i=0; i<size; i++)  BaryonField[IENum][i] *= rscale;
+  
+  //    remaining fields (all density-like)
+  int ifield;
+  rscale = DensityUnits;
+  for (ifield=0; ifield<DengoData.NumberOfFields; ifield++) {
+    for (i=0; i<size; i++) 
+      BaryonField[DengoData.field[ifield]][i] *= rscale;
+  }
+
+  //--------------
+  // call dengo
+  int ierr = 0;
+  
+  
+  //--------------
+  // rescale all relevant fields back to normalized units
+  //    internal energy
+  float rscale = LengthUnits*LengthUnits/TimeUnits/TimeUnits;
+  for (i=0; i<size; i++)  BaryonField[IENum][i] /= rscale;
+  
+  //    remaining fields (all density-like)
+  int ifield;
+  rscale = DensityUnits;
+  for (ifield=0; ifield<DengoData.NumberOfFields; ifield++) {
+    for (i=0; i<size; i++) 
+      BaryonField[DengoData.field[ifield]][i] /= rscale;
+  }
+
+
+  // check dengo error flag
+  if (ierr) {
+    fprintf(stdout, "GridLeftEdge = %"FSYM" %"FSYM" %"FSYM"\n",
+	    GridLeftEdge[0], GridLeftEdge[1], GridLeftEdge[2]);
+    fprintf(stdout, "GridRightEdge = %"FSYM" %"FSYM" %"FSYM"\n",
+	    GridRightEdge[0], GridRightEdge[1], GridRightEdge[2]);
+    fprintf(stdout, "GridDimension = %"ISYM" %"ISYM" %"ISYM"\n",
+	    GridDimension[0], GridDimension[1], GridDimension[2]);
+    ENZO_FAIL("Error in Dengo rate/cool solver!\n");
+  }
+
+
+
+  return SUCCESS;
+}