commit/enzo-dev: gbryan: Merged in ibutsky/enzo-dev (pull request #421)

[commits...@bitbucket.org]

1 new commit in enzo-dev:

https://bitbucket.org/enzo/enzo-dev/commits/4c537ac40bce/
Changeset: 4c537ac40bce
Branch: week-of-code
User: gbryan
Date: 2018-10-24 16:36:23+00:00
Summary: Merged in ibutsky/enzo-dev (pull request #421)

Fixing memory leak in magnetic supernova feedback

Approved-by: Britton Smith <britto...@gmail.com>
Approved-by: Brian O'Shea <os...@msu.edu>
Approved-by: Greg Bryan <gbr...@astro.columbia.edu>
Affected #: 28 files

diff -r bc76c0b9ebf0 -r 4c537ac40bce doc/manual/source/physics/star_particles.rst
--- a/doc/manual/source/physics/star_particles.rst
+++ b/doc/manual/source/physics/star_particles.rst
@@ -552,29 +552,36 @@
----------------------------
*Source: hydro_rk/SuperNovaSeedField.C*

-Select this method by setting ``UseSupernovaSeedFieldSourceTerms = 1``
-(Default = 0) and
-specifying the following parameters:
+Select this method by setting ``UseMagneticSupernovaFeedback = 1``
+(Default = 0) and specifying the parameters below. If
+``UseMagneticSupernovaFeedback == 2``, the parameters
+``MagneticSupernovaRadius`` and ``MagneticSupernovaDuration`` will
+be calculated to be the minimum allowed values (see below)
+at runtime based on the cell width and time step of the most-refined grid.

-``SupernovaSeedFieldTotalEnergy`` (in units of ergs) is the total amount
-of magnetic energy to be injected by a single supernova event. Defualt = 0.0.

-``SupernovaSeedFieldRadius`` (in units of parsecs) gives the scale over
+``MagneticSupernovaEnergy`` (in units of ergs) is the total amount
+of magnetic energy to be injected by a single supernova event. Defualt = 1e51
+ergs.
+
+``MagneticSupernovaRadius`` (in units of parsecs) gives the scale over
which to inject supernova energy. The injection mechanism normalizes the
spatial exponential decay of the injected supernova energy so that all of the
energy is contained within the specified radius. For this reason, the
-``SupernovaSeedFieldRadius`` should be at least 3 times the minimum cell width of
-the simulation. Default = 0.0.
+``MagneticSupernovaRadius`` should be at least 1.5 times the minimum cell width of
+the simulation (in pc). Default = 300 pc.

-``SupernovaSeedFieldDuration`` (in units of Myr) gives the duration of the
+``MagneticSupernovaDuration`` (in units of years) gives the duration of the
supernova magnetic energy injection. The injection mechanism is normalized so
-that all of the ``SupernovaSeedFieldTotalEnergy`` is injected over this
-time scale. In order to inject the correct amount of energy, ``SupernovaSeedFieldDuration`` should be set to at least 4
-times the minimum time step of the simulation. Default = 0.0.
+that all of the ``MagneticSupernovaEnergy`` is injected over this
+time scale. In order to inject the correct amount of energy, ``MagneticSupernovaDuration`` should be set to at least 5
+times the minimum time step of the simulation. Default = 50,000 years.


-The following applies to Methods 0 (Cen & Ostriker) and 1 (+
-stochastic star formation). The magnetic feedback method is described fully in `Butsky et al. (2017)
+The following applies to all star formation methods that produce a
+PARTICLE_TYPE_STAR object. Methods 0 (Cen & Ostriker) and 1 (+
+stochastic star formation) have been tested extensively.
+The magnetic feedback method is described fully in `Butsky et al. (2017)
<http://adsabs.harvard.edu/abs/2017ApJ...843..113B>`_.

When a star cluster particle reaches the end of its lifetime, we inject a

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/EnzoTiming.h
--- a/src/enzo/EnzoTiming.h
+++ b/src/enzo/EnzoTiming.h
@@ -29,8 +29,15 @@
#include <cstring>
#include <map>

-#define min(A,B) ((A) < (B) ? (A) : (B))
-#define max(A,B) ((A) > (B) ? (A) : (B))
+template <typename T>
+bool min(const T& A, const T& B) {
+ return A < B ? A : B;
+}
+
+template <typename T>
+bool max(const T& A, const T& B) {
+ return A > B ? A : B;
+}

double ReturnWallTime(void);
void Reduce_Times(double time, double *time_array);

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/EvolveLevel.C
--- a/src/enzo/EvolveLevel.C
+++ b/src/enzo/EvolveLevel.C
@@ -540,6 +540,10 @@

Grids[grid1]->GridData->CopyBaryonFieldToOldBaryonField();

+ // Find recently-supernova stars to add them the MagneticSupernovaList
+ if ((UseMagneticSupernovaFeedback) && (level == MaximumRefinementLevel))
+ Grids[grid1]->GridData->AddMagneticSupernovaeToList();
+
/* Call hydro solver and save fluxes around subgrids.
* HD_RK and MHD_RK are the 2nd order Runge-Kutta integrations, which
* require two steps (*_1stStep and *_2ndStep)
@@ -547,7 +551,7 @@
* All others (PPM, Zeus, MHD_Li/CT) are called from SolveHydroEquations
*/

-
+
if( HydroMethod != HD_RK && HydroMethod != MHD_RK ){
Grids[grid1]->GridData->SolveHydroEquations(LevelCycleCount[level],
NumberOfSubgrids[grid1], SubgridFluxesEstimate[grid1], level);
@@ -713,7 +717,10 @@
if (ComovingCoordinates)
Grids[grid1]->GridData->ComovingExpansionTerms();

- } // end loop over grids
+ if (UseMagneticSupernovaFeedback)
+ Grids[grid1]->GridData->MagneticSupernovaList.clear();
+
+ } // end loop over grids

/* Finalize (accretion, feedback, etc.) star particles */


diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/FreeRealMem.C
--- /dev/null
+++ b/src/enzo/FreeRealMem.C
@@ -0,0 +1,199 @@
+#ifdef TASKMAP
+
+#ifdef USE_LL
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <sys/types.h>
+#include "llapi.h"
+
+#include "macros_and_parameters.h"
+
+
+Eint64 FreeRealMem( char *node )
+{
+
+ LL_element *queryObject;
+ LL_element *job, *step;
+
+ Eint32 rc, obj_count, err_code, ok;
+
+ Eint64 free;
+ Eint64 freemem;
+ Eint64 def_mem;
+
+ def_mem = (Eint64)(12.0); /* minimum memory on DataStar */
+
+ Eint32 node_number;
+
+#ifdef SEABORG
+ sscanf(node, "s%5d", &node_number);
+ fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
+
+ if ( 701 <= node_number && node_number <= 713 ) {
+ def_mem = (Eint64)(56.0);
+ } else {
+
+ if ( 3701 <= node_number && node_number <= 5213 ) {
+ def_mem = (Eint64)(28.0);
+ } else {
+ def_mem = (Eint64)(14.0);
+ }
+
+ }
+#endif
+
+#ifdef DATASTAR
+ sscanf(node, "ds%3d", &node_number);
+ fprintf(stderr, "DataStar Node %s number %d\n", node, node_number);
+
+ if ( 4 <= node_number && node_number <= 11 ) {
+ def_mem = (Eint64)(120.0);
+ } else {
+
+ if ( 300 <= node_number && node_number <= 395 ) {
+ def_mem = (Eint64)(27.0);
+ } else {
+ def_mem = (Eint64)(12.0);
+ }
+
+ }
+#endif
+
+ def_mem = def_mem * 1024.0;
+ fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
+
+ ok = 0;
+
+ /* Initialize the query */
+
+ queryObject = ll_query(MACHINES);
+
+ if (!queryObject) {
+ fprintf(stderr, "Query MACHINES: ll_query() returns NULL.\n");
+ freemem = def_mem;
+ return(freemem);
+ }
+
+ char **nodelist;
+ nodelist = (char **) malloc(2*sizeof(char *));
+ char nodename[6];
+ strcpy(nodename, node);
+
+ nodelist[0] = nodename;
+ nodelist[1] = NULL;
+
+ rc = ll_set_request(queryObject, QUERY_HOST, nodelist, ALL_DATA);
+
+ if (rc) {
+ fprintf(stderr, "Query MACHINES: ll_set_request() return code is non-zero.\n");
+ freemem = def_mem;
+ return(freemem);
+ }
+
+ job = ll_get_objs(queryObject, LL_CM, NULL, &obj_count, &err_code);
+
+ /* Process the job objects */
+
+ rc = ll_get_data(job, LL_MachineFreeRealMemory64, &free);
+
+ if (rc) {
+ fprintf(stderr, "Get data: ll_get_data() return code is non-zero: %d\n", rc);
+ freemem = def_mem;
+ return(freemem);
+ }
+
+/*
+ fprintf(stderr, "Free: %lld MBytes\n", free);
+*/
+
+ ll_free_objs(queryObject);
+ ll_deallocate(queryObject);
+
+ freemem = (Eint64)(free);
+
+ return(freemem);
+}
+
+#else
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "macros_and_parameters.h"
+
+Eint64 FreeRealMem( char *node )
+{
+
+ Eint64 freemem;
+ Eint64 def_mem;
+
+ def_mem = (Eint64)(12.5); /* minimum memory of DataStar */
+
+ Eint32 node_number;
+
+
+#ifdef SEABORG
+ sscanf(node, "s%5d", &node_number);
+ fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
+
+ if ( 701 <= node_number && node_number <= 713 ) {
+ def_mem = (Eint64)(56.0);
+ } else {
+
+ if ( 3701 <= node_number && node_number <= 5213 ) {
+ def_mem = (Eint64)(28.0);
+ } else {
+ def_mem = (Eint64)(14.0);
+ }
+
+ }
+#endif
+
+#ifdef DATASTAR
+ sscanf(node, "ds%3d", &node_number);
+ fprintf(stderr, "DataStarNode %s number %d\n", node, node_number);
+
+ if ( 4 <= node_number && node_number <= 11 ) {
+ def_mem = (Eint64)(120.0);
+ } else {
+
+ if ( 300 <= node_number && node_number <= 395 ) {
+ def_mem = (Eint64)(27.0);
+ } else {
+ def_mem = (Eint64)(12.0);
+ }
+
+ }
+#endif
+
+ def_mem = def_mem * 1024.0;
+ fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
+
+ freemem = def_mem;
+
+ return(freemem);
+}
+
+#endif
+
+#else
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "macros_and_parameters.h"
+
+Eint64 FreeRealMem( char *node )
+{
+
+ Eint64 freemem;
+
+ fprintf(stderr, "Error in FreeRealMem - the TASKMAP macro is not defined! Default to 12.0 GBytes.\n");
+
+ freemem = (Eint64)(12.0*1024.0); /* minimum memory on DataStar */
+
+ return(freemem);
+}
+
+#endif

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/FreeRealMem.c
--- a/src/enzo/FreeRealMem.c
+++ /dev/null
@@ -1,199 +0,0 @@
-#ifdef TASKMAP
-
-#ifdef USE_LL
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <sys/types.h>
-#include "llapi.h"
-
-#include "macros_and_parameters.h"
-
-
-Eint64 FreeRealMem( char *node )
-{
-
- LL_element *queryObject;
- LL_element *job, *step;
-
- Eint32 rc, obj_count, err_code, ok;
-
- Eint64 free;
- Eint64 freemem;
- Eint64 def_mem;
-
- def_mem = (Eint64)(12.0); /* minimum memory on DataStar */
-
- Eint32 node_number;
-
-#ifdef SEABORG
- sscanf(node, "s%5d", &node_number);
- fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
-
- if ( 701 <= node_number && node_number <= 713 ) {
- def_mem = (Eint64)(56.0);
- } else {
-
- if ( 3701 <= node_number && node_number <= 5213 ) {
- def_mem = (Eint64)(28.0);
- } else {
- def_mem = (Eint64)(14.0);
- }
-
- }
-#endif
-
-#ifdef DATASTAR
- sscanf(node, "ds%3d", &node_number);
- fprintf(stderr, "DataStar Node %s number %d\n", node, node_number);
-
- if ( 4 <= node_number && node_number <= 11 ) {
- def_mem = (Eint64)(120.0);
- } else {
-
- if ( 300 <= node_number && node_number <= 395 ) {
- def_mem = (Eint64)(27.0);
- } else {
- def_mem = (Eint64)(12.0);
- }
-
- }
-#endif
-
- def_mem = def_mem * 1024.0;
- fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
-
- ok = 0;
-
- /* Initialize the query */
-
- queryObject = ll_query(MACHINES);
-
- if (!queryObject) {
- fprintf(stderr, "Query MACHINES: ll_query() returns NULL.\n");
- freemem = def_mem;
- return(freemem);
- }
-
- char **nodelist;
- nodelist = (char **) malloc(2*sizeof(char *));
- char nodename[6];
- strcpy(nodename, node);
-
- nodelist[0] = nodename;
- nodelist[1] = NULL;
-
- rc = ll_set_request(queryObject, QUERY_HOST, nodelist, ALL_DATA);
-
- if (rc) {
- fprintf(stderr, "Query MACHINES: ll_set_request() return code is non-zero.\n");
- freemem = def_mem;
- return(freemem);
- }
-
- job = ll_get_objs(queryObject, LL_CM, NULL, &obj_count, &err_code);
-
- /* Process the job objects */
-
- rc = ll_get_data(job, LL_MachineFreeRealMemory64, &free);
-
- if (rc) {
- fprintf(stderr, "Get data: ll_get_data() return code is non-zero: %d\n", rc);
- freemem = def_mem;
- return(freemem);
- }
-
-/*
- fprintf(stderr, "Free: %lld MBytes\n", free);
-*/
-
- ll_free_objs(queryObject);
- ll_deallocate(queryObject);
-
- freemem = (Eint64)(free);
-
- return(freemem);
-}
-
-#else
-
-#include <stdlib.h>
-#include <stdio.h>
-
-#include "macros_and_parameters.h"
-
-Eint64 FreeRealMem( char *node )
-{
-
- Eint64 freemem;
- Eint64 def_mem;
-
- def_mem = (Eint64)(12.5); /* minimum memory of DataStar */
-
- Eint32 node_number;
-
-
-#ifdef SEABORG
- sscanf(node, "s%5d", &node_number);
- fprintf(stderr, "Seaborg Node %s number %d\n", node, node_number);
-
- if ( 701 <= node_number && node_number <= 713 ) {
- def_mem = (Eint64)(56.0);
- } else {
-
- if ( 3701 <= node_number && node_number <= 5213 ) {
- def_mem = (Eint64)(28.0);
- } else {
- def_mem = (Eint64)(14.0);
- }
-
- }
-#endif
-
-#ifdef DATASTAR
- sscanf(node, "ds%3d", &node_number);
- fprintf(stderr, "DataStarNode %s number %d\n", node, node_number);
-
- if ( 4 <= node_number && node_number <= 11 ) {
- def_mem = (Eint64)(120.0);
- } else {
-
- if ( 300 <= node_number && node_number <= 395 ) {
- def_mem = (Eint64)(27.0);
- } else {
- def_mem = (Eint64)(12.0);
- }
-
- }
-#endif
-
- def_mem = def_mem * 1024.0;
- fprintf(stderr,"Default node memory on node %s [%d] is %"ISYM" MBytes\n", node, node_number, def_mem);
-
- freemem = def_mem;
-
- return(freemem);
-}
-
-#endif
-
-#else
-
-#include <stdlib.h>
-#include <stdio.h>
-
-#include "macros_and_parameters.h"
-
-Eint64 FreeRealMem( char *node )
-{
-
- Eint64 freemem;
-
- fprintf(stderr, "Error in FreeRealMem - the TASKMAP macro is not defined! Default to 12.0 GBytes.\n");
-
- freemem = (Eint64)(12.0*1024.0); /* minimum memory on DataStar */
-
- return(freemem);
-}
-
-#endif

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Grid.h
--- a/src/enzo/Grid.h
+++ b/src/enzo/Grid.h
@@ -11,6 +11,7 @@
************************************************************************/
#ifndef GRID_DEFINED__
#define GRID_DEFINED__
+#include <vector>
#include "ProtoSubgrid.h"
#include "ListOfParticles.h"
#include "region.h"
@@ -20,7 +21,6 @@
#include "Star.h"
#include "FOF_allvars.h"
#include "MemoryPool.h"
-#include "list.h"
#include "hydro_rk/SuperNova.h"
#ifdef ECUDA
#include "hydro_rk/CudaMHD.h"
@@ -3011,11 +3011,10 @@
int MHDCT_ConvertEnergyToConservedS();
int MHDCT_ConvertEnergyToSpecificS();

- //List of SuperNova objects that each grid needs to keep track of
-
- List<SuperNova> SuperNovaList;
-
-
+ // used if UseMagneticSupernovaFeedback is turned on
+ int AddMagneticSupernovaeToList();
+ //List of SuperNova objects that each grid needs to keep track of
+ std::vector<SuperNova> MagneticSupernovaList;
};

// inline int grid::ReadRandomForcingFields (FILE *main_file_pointer);

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Grid_AddFeedbackSphere.C
--- a/src/enzo/Grid_AddFeedbackSphere.C
+++ b/src/enzo/Grid_AddFeedbackSphere.C
@@ -43,8 +43,6 @@

int FindField(int field, int farray[], int numfields);

-void mt_init(unsigned_int seed);
-
int grid::AddFeedbackSphere(Star *cstar, int level, float radius, float DensityUnits,
float LengthUnits, float VelocityUnits,
float TemperatureUnits, float TimeUnits, double EjectaDensity,
@@ -990,48 +988,7 @@
/* Now it's done, unmark. */

//cstar->FeedbackFlag = NO_FEEDBACK;
-
-
- // Create a randomly-oriented SuperNova object and add it to the SuperNova Grid list
- if (cstar->FeedbackFlag == SUPERNOVA_SEEDFIELD) {
- if(UseSupernovaSeedFieldSourceTerms){
-
- SuperNova *P = new SuperNova();
- int mt_seed = 100;
- mt_init(mt_seed);
- float random_u = (float)mt_seed/100.0; // random variable from 0 to 1
- float random_v = (float)mt_seed/100.0;
- float random_phi = 2*M_PI*random_u; // 0 to 2pi
- float random_theta = acos(2*random_v-1); // 0 to pi
- // Setting up randomly oriented magnetic feedback of supernova
- float phi_x = sin(random_theta)*cos(random_phi);
- float phi_y = sin(random_theta)*sin(random_phi);
- float phi_z = cos(random_theta);
-
- // Birthtime of supernova is at the end of a star particle's life
- float sn_birthtime = cstar->BirthTime + cstar->LifeTime;
-
- // Convert units to system units
- // Converting time from years to seconds, then internal units
- float sn_duration = SupernovaSeedFieldDuration * 3.1556952e7 / TimeUnits;
- // Converting radius from parsecs to cm, then internal units
- float sn_radius = SupernovaSeedFieldRadius * 3.0856775714e18 / LengthUnits;
- // Converting energy from ergs to internal units
- float MassUnits = DensityUnits * POW(LengthUnits, 3);
- float sn_energy = SupernovaSeedFieldEnergy / (MassUnits*VelocityUnits*VelocityUnits);
-
- // Creates a supernova with magnetic feedback set by user-defined parameters and
- // adds it to supernova list
- if((Time > sn_birthtime) && (Time < sn_birthtime + SupernovaSeedFieldDuration)){
- P->setValues(phi_x, phi_y, phi_z, cstar->pos[0], cstar->pos[1], cstar->pos[2],
- sn_radius, sn_birthtime, sn_duration, sn_energy);
-
- this->SuperNovaList.append(P);
- printf("Added supernova to list \n");
- }
- }
- }
-
+

return SUCCESS;


diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Grid_AddMagneticSupernovaeToList.C
--- /dev/null
+++ b/src/enzo/Grid_AddMagneticSupernovaeToList.C
@@ -0,0 +1,112 @@
+/***********************************************************************
+/
+/ FIND DEAD STAR PARTICLES AND ADD THEM TO MAGNETIC SUPERNOVA LIST
+/
+/ written by: Iryna Butsky
+/ date: June 2018
+/
+/ PURPOSE: To apply magnetic supernova feedback, we need to first add
+/ all stars that have recently gone supernova to the global list.
+/
+************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <math.h>
+
+#include "ErrorExceptions.h"
+#include "macros_and_parameters.h"
+#include "typedefs.h"
+#include "global_data.h"
+#include "phys_constants.h"
+#include "Fluxes.h"
+#include "GridList.h"
+#include "ExternalBoundary.h"
+#include "TopGridData.h"
+#include "Grid.h"
+#include "hydro_rk/SuperNova.h"
+
+int GetUnits(float *DensityUnits, float *LengthUnits,
+ float *TemperatureUnits, float *TimeUnits,
+ float *VelocityUnits, double *MassUnits, FLOAT Time);
+
+void mt_init(unsigned_int seed);
+unsigned_long_int mt_random();
+
+int grid::AddMagneticSupernovaeToList()
+{
+
+ if (ProcessorNumber != MyProcessorNumber) {
+ return SUCCESS;
+ }
+
+ float random_u, random_v, random_phi, random_theta, phi_x, phi_y, phi_z;
+ float sn_birthtime, sn_duration, sn_radius, sn_energy, star_birthtime, star_lifetime;
+
+ float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, VelocityUnits, MassUnits, EnergyUnits;
+ if (GetUnits(&DensityUnits, &LengthUnits,&TemperatureUnits, &TimeUnits,
+ &VelocityUnits, &MassUnits, Time) == FAIL){
+ fprintf(stderr, "Error in GetUnits.\n");
+ return FAIL;
+ }
+
+
+ // if UseMagneticSupernovaFeedback > 1, then we set the magnetic feedback radius and duration
+ // below based on the resolution of the grid at the highest refinement level
+ if (UseMagneticSupernovaFeedback > 1){
+ sn_duration = 5.0 * dtFixed;
+ MagneticSupernovaDuration = sn_duration * TimeUnits / SecondsPerYear;
+ sn_radius = 1.5 * this->CellWidth[0][0];
+ MagneticSupernovaRadius = sn_radius * LengthUnits / pc;
+ }
+ else {
+ // Converting time from years to seconds, then internal units
+ sn_duration = MagneticSupernovaDuration * SecondsPerYear / TimeUnits;
+ if(sn_duration < 5.0 * dtFixed) {
+ printf("WARNING: Magnetic supernova feedback duration is less than the recommended minimum of 5 x dtFixed\n");
+ printf("Current dtFixed = %e years \n", dtFixed * TimeUnits / SecondsPerYear);
+ }
+ // Converting radius from parsecs to cm, then internal units
+ sn_radius = MagneticSupernovaRadius * pc / LengthUnits;
+ if(sn_radius < 1.5 * this->CellWidth[0][0]){
+ printf("WARNING: Magnetic supernova feedback radius is less than the recommended minimum of 1.5 x CellWidth\n");
+ printf("Current CellWidth = %e pc \n", this->CellWidth[0][0] * LengthUnits / pc);
+ }
+
+ }
+ // Converting energy from ergs to internal units
+ sn_energy = MagneticSupernovaEnergy / (MassUnits*VelocityUnits*VelocityUnits);
+
+ for (int star_id = 0; star_id < this->NumberOfParticles; star_id++){
+ if (this->ParticleType[star_id] == PARTICLE_TYPE_STAR){
+
+ star_birthtime = this->ParticleAttribute[0][star_id];
+ star_lifetime = this->ParticleAttribute[1][star_id];
+ sn_birthtime = star_birthtime + star_lifetime;
+
+ // Creates a supernova with magnetic feedback set by user-defined parameters and
+ // adds it to supernova list
+ if((this->Time > sn_birthtime) && (this->Time < sn_birthtime + sn_duration)){
+ SuperNova P = SuperNova();
+ mt_init((unsigned_int) this->ParticleNumber[star_id]);
+
+ random_u = (float)(mt_random()%32768)/32768.0; // random variable from 0 to 1
+ random_v = (float)(mt_random()%32768)/32768.0;
+ random_phi = 2*M_PI*random_u; // 0 to 2pi
+ random_theta = acos(2*random_v-1); // 0 to pi
+
+ // Setting up randomly oriented magnetic feedback of supernova
+ phi_x = sin(random_theta)*cos(random_phi);
+ phi_y = sin(random_theta)*sin(random_phi);
+ phi_z = cos(random_theta);
+
+ P.setValues(phi_x, phi_y, phi_z,
+ ParticlePosition[0][star_id], ParticlePosition[1][star_id], ParticlePosition[2][star_id],
+ sn_radius, sn_birthtime, sn_duration, sn_energy);
+ this->MagneticSupernovaList.push_back(P);
+ }
+ }
+ }
+ return SUCCESS;
+}

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Grid_FindAllStarParticles.C
--- a/src/enzo/Grid_FindAllStarParticles.C
+++ b/src/enzo/Grid_FindAllStarParticles.C
@@ -50,7 +50,6 @@
//StarType = abs(ParticleType[i]);
if (ParticleType[i] == PARTICLE_TYPE_SINGLE_STAR ||
ParticleType[i] == PARTICLE_TYPE_BLACK_HOLE ||
- (ParticleType[i] == PARTICLE_TYPE_STAR && UseSupernovaSeedFieldSourceTerms) ||
ParticleType[i] == PARTICLE_TYPE_CLUSTER ||
ParticleType[i] == PARTICLE_TYPE_COLOR_STAR ||
ParticleType[i] == PARTICLE_TYPE_MBH ||
@@ -58,7 +57,7 @@
(StarParticleRadiativeFeedback == TRUE &&
ParticleType[i] == PARTICLE_TYPE_STAR)) {

- if (ParticleType[i] == PARTICLE_TYPE_STAR)
+ if (ParticleType[i] == PARTICLE_TYPE_STAR && UseMagneticSupernovaFeedback == 0)
LifetimeFactor = 12.0;
else
LifetimeFactor = 1.0;

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Grid_FindNewStarParticles.C
--- a/src/enzo/Grid_FindNewStarParticles.C
+++ b/src/enzo/Grid_FindNewStarParticles.C
@@ -39,9 +39,6 @@

for (i = 0; i < NumberOfParticles; i++)
if (ParticleType[i] == -PARTICLE_TYPE_SINGLE_STAR ||
- (ParticleType[i] == -PARTICLE_TYPE_STAR && UseSupernovaSeedFieldSourceTerms &&
- (this->Time >= ParticleAttribute[0][i] &&
- this->Time <= ParticleAttribute[0][i]+ParticleAttribute[1][i])) ||
ParticleType[i] == -PARTICLE_TYPE_BLACK_HOLE ||
ParticleType[i] == -PARTICLE_TYPE_CLUSTER ||
ParticleType[i] == -PARTICLE_TYPE_COLOR_STAR ||

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Grid_constructor.C
--- a/src/enzo/Grid_constructor.C
+++ b/src/enzo/Grid_constructor.C
@@ -200,12 +200,11 @@
MagneticDims[field][dim] = -100;
}}

-
/* For once-per-rootgrid-timestep star formation, the following flag
determines whether SF is about to occur or not. It's currently
(April 2012) only implemented for H2REG_STAR and completely
ignored for all other star makers. */
MakeStars = 0;
- if (UseSupernovaSeedFieldSourceTerms == 1) List<SuperNova> SuperNovaList;

+ if (UseMagneticSupernovaFeedback) std::vector<SuperNova> MagneticSupernovaList;
}

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Make.config.objects
--- a/src/enzo/Make.config.objects
+++ b/src/enzo/Make.config.objects
@@ -307,6 +307,7 @@
Grid_AddFeedbackSphere.o \
Grid_AddFieldMassToMassFlaggingField.o \
Grid_AddFields.o \
+ Grid_AddMagneticSupernovaeToList.o \
Grid_AddOneParticleFromList.o \
Grid_AddOverlappingParticleMassField.o \
Grid_AddParticlesFromList.o \

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/ReadParameterFile.C
--- a/src/enzo/ReadParameterFile.C
+++ b/src/enzo/ReadParameterFile.C
@@ -1273,10 +1273,10 @@
ret += sscanf(line, "GasDragCoefficient = %"GSYM, &GasDragCoefficient);

// Parameters for magnetic feedback from supernovae
- ret += sscanf(line, "UseSupernovaSeedFieldSourceTerms = %"ISYM, &UseSupernovaSeedFieldSourceTerms);
- ret += sscanf(line,"SupernovaSeedFieldRadius = %"FSYM, &SupernovaSeedFieldRadius);
- ret += sscanf(line,"SupernovaSeedFieldEnergy = %"FSYM, &SupernovaSeedFieldEnergy);
- ret += sscanf(line,"SupernovaSeedFieldDuration = %"FSYM, &SupernovaSeedFieldDuration);
+ ret += sscanf(line, "UseMagneticSupernovaFeedback = %"ISYM, &UseMagneticSupernovaFeedback);
+ ret += sscanf(line,"MagneticSupernovaRadius = %"FSYM, &MagneticSupernovaRadius);
+ ret += sscanf(line,"MagneticSupernovaEnergy = %"FSYM, &MagneticSupernovaEnergy);
+ ret += sscanf(line,"MagneticSupernovaDuration = %"FSYM, &MagneticSupernovaDuration);

/* If the dummy char space was used, then make another. */


diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/SetDefaultGlobalValues.C
--- a/src/enzo/SetDefaultGlobalValues.C
+++ b/src/enzo/SetDefaultGlobalValues.C
@@ -978,12 +978,12 @@
UseGasDrag = 0;
GasDragCoefficient = 0.;

- /* Supernova magnetic seed field */
+ /* Supernova magnetic feedback */
/* Default == 0 -> no magnetic field contribution */
- UseSupernovaSeedFieldSourceTerms = 0;
- SupernovaSeedFieldRadius = 0.0;
- SupernovaSeedFieldDuration = 0.0;
- SupernovaSeedFieldEnergy = 0.0;
+ UseMagneticSupernovaFeedback = 0; // 0 = off; 1+ = on; 2 = Radius and duration calculated during runtime
+ MagneticSupernovaRadius = 300.0; // Total injection radius of magnetic field in parsecs
+ MagneticSupernovaDuration = 5e4 ; // Total duration of magnetic feedback in years
+ MagneticSupernovaEnergy = 1.0e51; // Total energy (ergs) injected per star particle (supernova)

return SUCCESS;
}

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/StarParticleAddFeedback.C
--- a/src/enzo/StarParticleAddFeedback.C
+++ b/src/enzo/StarParticleAddFeedback.C
@@ -88,10 +88,7 @@
&TimeUnits, &VelocityUnits, Time);

count = 0;
- // clear list of Supernovae at each timestep to avoid adding duplicates in Grid_AddFeedbackSphere
- if(UseSupernovaSeedFieldSourceTerms){
- LevelArray[level]->GridData->SuperNovaList.clear();
- }
+
for (cstar = AllStars; cstar; cstar = cstar->NextStar, count++) {

AddedFeedback[count] = false;

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Star_Accrete.C
--- a/src/enzo/Star_Accrete.C
+++ b/src/enzo/Star_Accrete.C
@@ -29,15 +29,10 @@

int Star::Accrete(void)
{
- if (UseSupernovaSeedFieldSourceTerms == 1){
- if (this->CurrentGrid == NULL ||(this->naccretions == 0))
- return SUCCESS;
- }
- else {
- if (this->CurrentGrid == NULL ||
- (this->naccretions == 0 && fabs(this->DeltaMass) < tiny_number))
- return SUCCESS;
- }
+ if (this->CurrentGrid == NULL ||
+ (this->naccretions == 0 && fabs(this->DeltaMass) < tiny_number))
+ return SUCCESS;
+

const double Msun = 1.989e33, yr = 3.1557e7;
int dim, i, n, count;

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Star_HitEndpoint.C
--- a/src/enzo/Star_HitEndpoint.C
+++ b/src/enzo/Star_HitEndpoint.C
@@ -101,9 +101,6 @@
case PopIII_CF:
break;

- case NormalStar:
- break;
-
} // ENDSWITCH

return result;

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/Star_SetFeedbackFlag.C
--- a/src/enzo/Star_SetFeedbackFlag.C
+++ b/src/enzo/Star_SetFeedbackFlag.C
@@ -135,12 +135,7 @@

//this->type = abs_type;

- case PARTICLE_TYPE_STAR:
- if(UseSupernovaSeedFieldSourceTerms)
- this->FeedbackFlag = SUPERNOVA_SEEDFIELD;
- break;
-}
-
+ }

return SUCCESS;
}

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/WriteParameterFile.C
--- a/src/enzo/WriteParameterFile.C
+++ b/src/enzo/WriteParameterFile.C
@@ -1200,10 +1200,10 @@
fprintf(fptr, "CorrectParentBoundaryFlux = %d\n", CorrectParentBoundaryFlux);

/* Supernova magnetic seed field */
- fprintf(fptr, "UseSupernovaSeedFieldSourceTerms = %"ISYM"\n", UseSupernovaSeedFieldSourceTerms);
- fprintf(fptr, "SupernovaSeedFieldRadius = %"GSYM"\n", SupernovaSeedFieldRadius);
- fprintf(fptr,"SupernovaSeedFieldEnergy = %"GSYM"\n",SupernovaSeedFieldEnergy);
- fprintf(fptr,"SupernovaSeedFieldDuration = %"GSYM"\n",SupernovaSeedFieldDuration);
+ fprintf(fptr, "UseMagneticSupernovaFeedback = %"ISYM"\n", UseMagneticSupernovaFeedback);
+ fprintf(fptr, "MagneticSupernovaRadius = %"GSYM"\n", MagneticSupernovaRadius);
+ fprintf(fptr,"MagneticSupernovaEnergy = %"GSYM"\n",MagneticSupernovaEnergy);
+ fprintf(fptr,"MagneticSupernovaDuration = %"GSYM"\n",MagneticSupernovaDuration);

/* Output current time */
time_t ID;

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/global_data.h
--- a/src/enzo/global_data.h
+++ b/src/enzo/global_data.h
@@ -1109,9 +1109,9 @@
EXTERN float GalaxySimulationPreWindVelocity[MAX_DIMENSION];

/* Supernova magnetic seed field */
-EXTERN int UseSupernovaSeedFieldSourceTerms;
-EXTERN float SupernovaSeedFieldRadius;
-EXTERN float SupernovaSeedFieldDuration;
-EXTERN float SupernovaSeedFieldEnergy;
+EXTERN int UseMagneticSupernovaFeedback;
+EXTERN float MagneticSupernovaRadius;
+EXTERN float MagneticSupernovaDuration;
+EXTERN float MagneticSupernovaEnergy;

#endif

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/hydro_rk/Grid_MHDSourceTerms.C
--- a/src/enzo/hydro_rk/Grid_MHDSourceTerms.C
+++ b/src/enzo/hydro_rk/Grid_MHDSourceTerms.C
@@ -19,13 +19,14 @@
#include "macros_and_parameters.h"
#include "typedefs.h"
#include "global_data.h"
-#include "list.h"
+#include "phys_constants.h"
#include "Fluxes.h"
#include "GridList.h"
#include "ExternalBoundary.h"
#include "TopGridData.h"
#include "Grid.h"
#include "EOS.h"
+#include "hydro_rk/SuperNova.h"

int GetUnits(float *DensityUnits, float *LengthUnits,
float *TemperatureUnits, float *TimeUnits,
@@ -33,6 +34,8 @@
int CosmologyComputeExpansionFactor(FLOAT time, FLOAT *a, FLOAT *dadt);
int FindField(int field, int farray[], int numfields);

+void mt_init(unsigned_int seed);
+unsigned_long_int mt_random();

int grid::MHDSourceTerms(float **dU)
{
@@ -357,49 +360,73 @@
}
}

-if((UseSupernovaSeedFieldSourceTerms == 1)) {
+ if(UseMagneticSupernovaFeedback) {
+
+ int n, active_x, active_y, center_i, center_j, center_k, num_sn_cells_x, num_sn_cells_y, num_sn_cells_z;
+ snsf_source_terms S;
+ float dx, dy, dz, dist_to_sn, magnetic_energy_density;
+ float DensityUnits, LengthUnits, TemperatureUnits, TimeUnits, VelocityUnits;
+
+ if (GetUnits(&DensityUnits, &LengthUnits,&TemperatureUnits, &TimeUnits,
+ &VelocityUnits, Time) == FAIL){
+ fprintf(stderr, "Error in GetUnits.\n");
+ return FAIL;
+ }
+
+ // Converting radius from parsecs to cm, then internal units
+ float sn_radius = MagneticSupernovaRadius * pc / LengthUnits;

- int n = 0, igrid;
- int iden=FindField(Density, FieldType, NumberOfBaryonFields);
- snsf_source_terms S;
- List<SuperNova>::Iterator *P = this->SuperNovaList.begin();
- FLOAT cell_center[3];
- FLOAT dx, dy, dz, dist_to_sn;
- int temp =1;
- int entered = 0;
+ active_x = GridDimension[0] - 2*NumberOfGhostZones;
+ active_y = GridDimension[1] - 2*NumberOfGhostZones;

- for (int k = GridStartIndex[2]; k <= GridEndIndex[2]; k++) {
- for (int j = GridStartIndex[1]; j <= GridEndIndex[1]; j++) {
- for (int i = GridStartIndex[0]; i <= GridEndIndex[0]; i++, n++) {
- igrid = i+(j+k*GridDimension[1])*GridDimension[0];
+ // the number of additional cells that receive magnetic supernova feedback.
+ // assuming the supernova is in the center of one cell
+ num_sn_cells_x = (int) ((sn_radius - 0.5*CellWidth[0][0])/ CellWidth[0][0]);
+ num_sn_cells_y = (int) ((sn_radius - 0.5*CellWidth[1][0])/ CellWidth[1][0]);
+ num_sn_cells_z = (int) ((sn_radius - 0.5*CellWidth[2][0])/ CellWidth[2][0]);
+
+ for (std::vector<SuperNova>::iterator current_sn = this->MagneticSupernovaList.begin();
+ current_sn != this->MagneticSupernovaList.end(); current_sn++){

- cell_center[0] = CellLeftEdge[0][i] + 0.5*CellWidth[0][i];
- cell_center[1] = CellLeftEdge[1][j] + 0.5*CellWidth[1][j];
- cell_center[2] = CellLeftEdge[2][k] + 0.5*CellWidth[2][k];
+ // find index of the cell nearest to the supernova center
+ // assuming that supernova in the center of that cell
+ center_i = (int)((current_sn->getPosition()[0] - GridLeftEdge[0]) / CellWidth[0][0]);
+ center_j = (int)((current_sn->getPosition()[1] - GridLeftEdge[1]) / CellWidth[1][0]);
+ center_k = (int)((current_sn->getPosition()[2] - GridLeftEdge[2]) / CellWidth[2][0]);

- while(P != this->SuperNovaList.end()){
- dx = P->get()->getPosition()[0] - cell_center[0];
- dy = P->get()->getPosition()[1] - cell_center[1];
- dz = P->get()->getPosition()[2] - cell_center[2];
-
- dist_to_sn = sqrt(dx*dx + dy*dy + dz*dz);
- if (dist_to_sn < 1.1*SupernovaSeedFieldRadius){
- S = P->get()->getSourceTerms(dx, dy, dz, Time);
- double rho = BaryonField[DensNum][igrid];
+ for(int k = center_k - num_sn_cells_z; k <= center_k + num_sn_cells_z; k++){
+ for(int j = center_j - num_sn_cells_y; j <= center_j + num_sn_cells_y; j++){
+ for(int i = center_i - num_sn_cells_x; i <= center_i + num_sn_cells_x; i++){
+
+ // only add magnetic feedback on the active grid cells
+ if ((k >= GridStartIndex[2]) && (k <= GridEndIndex[2]) &&
+ (j >= GridStartIndex[1]) && (j <= GridEndIndex[1]) &&
+ (i >= GridStartIndex[0]) && (i <= GridEndIndex[0])){
+
+ dx = CellWidth[0][0] * (float)(i-center_i);
+ dy = CellWidth[1][0] * (float)(j-center_j);
+ dz = CellWidth[2][0] * (float)(k-center_k);
+
+ dist_to_sn = sqrt(dx*dx + dy*dy + dz*dz);
+ S = current_sn->getSourceTerms(dx, dy, dz, Time);
+
+ // solving for index n
+ // analogous to how igrid is calculated, but taking into acount Ghost Zones
+ n = (i - GridStartIndex[0])+((j-GridStartIndex[1]) + (k-GridStartIndex[2])*active_y) * active_x;
+
+ dU[iBx][n] += S.dbx*dtFixed;
+ dU[iBy][n] += S.dby*dtFixed;
+ dU[iBz][n] += S.dbz*dtFixed;

- dU[iBx][n] += S.dbx*dtFixed;
- dU[iBy][n] += S.dby*dtFixed;
- dU[iBz][n] += S.dbz*dtFixed;
- dU[iEtot][n] += S.dUtot*dtFixed;
+ dU[iEtot][n] += S.dUtot * dtFixed;
+
+ }

- }
- P = P->next();
- }// End of SuperNovaList iteration
- } // End of k for-loop
- } // End of j for-loop
- } // End of i for-loop
-
-} // End of UseSuperNovaSeedFieldSourceTerms scope
+ } // End of k for-loop
+ } // End of j for-loop
+ } // End of i for-loop
+ } // End of MagneticSupernovaList loop
+ } // End of UseMagneticSupernovaFeedback scope




diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/hydro_rk/LLF_Zero_MHD.C
--- a/src/enzo/hydro_rk/LLF_Zero_MHD.C
+++ b/src/enzo/hydro_rk/LLF_Zero_MHD.C
@@ -34,7 +34,7 @@

// compute priml and primr
int iprim;
- for (int field = 0; field < NEQ_MHD; field++) {
+ for (int field = 0; field < NEQ_MHD-DualEnergyFormalism; field++) {
for (int i = 0; i < ActiveSize+1; i++) {
iprim = i + NumberOfGhostZones - 1;
priml[field][i] = prim[field][iprim];
@@ -53,9 +53,9 @@
iprim = i+NumberOfGhostZones-1;
for (int field = 0; field < NSpecies; field++) {
if (FluxLine[iD][i] >= 0) {
- species[field][i] = prim[field+5][iprim ];
+ species[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism][iprim ];
} else {
- species[field][i] = prim[field+5][iprim+1];
+ species[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism][iprim+1];
}
}
sum = 0;
@@ -76,9 +76,9 @@
iprim = i+NumberOfGhostZones-1;
for (int field = 0; field < NColor; field++) {
if (FluxLine[iD][i] >= 0) {
- colors[field][i] = prim[field+5+NSpecies][iprim ];
+ colors[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism+NSpecies][iprim ];
} else {
- colors[field][i] = prim[field+5+NSpecies][iprim+1];
+ colors[field][i] = prim[field+NEQ_MHD-DualEnergyFormalism+NSpecies][iprim+1];
}
}
for (int field = 0; field < NColor; field++) {

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/hydro_rk/SuperNovaSeedField.C
--- a/src/enzo/hydro_rk/SuperNovaSeedField.C
+++ b/src/enzo/hydro_rk/SuperNovaSeedField.C
@@ -127,14 +127,14 @@

snsf_source_terms SuperNova::getSourceTerms(double dx, double dy, double dz, double enzoTime){
// ------------------------------------------------------------------- //
- // Inputs: (1) distances to current cell, dx, dy, dz //
+ // Inputs: (1) distances to current cell, dx, dy, dz //
// (2) current time of the simulation //
// //
- // Output: (1) A struct of source terms corresponding to the given //
- // supernova //
- // (i) x,y,z components of time derivative of B-field //
- // (ii) time derivative of magnetic energy density //
- // (iii) time derivative of total energy density //
+ // Output: (1) A struct of source terms corresponding to the given //
+ // supernova //
+ // (i) x,y,z components of time derivative of B-field //
+ // (ii) time derivative of magnetic energy density //
+ // (iii) time derivative of total energy density //
// ------------------------------------------------------------------- //

FLOAT phi[3], zhat_cell[3];
@@ -180,6 +180,9 @@

FLOAT r_scale = (r_cyl/characteristicLength)*\
exp(- r_s*r_s /(characteristicLength*characteristicLength));
+
+ FLOAT db_r_scale = (r_cyl/characteristicLength)*\
+ exp(- r_s*r_s /(2.0*characteristicLength*characteristicLength));
FLOAT t_exp = exp(-t_s/characteristicTime);

FLOAT db_t_exp = t_exp / characteristicTime;

diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/macros_and_parameters.h
--- a/src/enzo/macros_and_parameters.h
+++ b/src/enzo/macros_and_parameters.h
@@ -1,5 +1,12 @@
#ifndef __macros_and_parameters_h_
#define __macros_and_parameters_h_
+#include "string.h"
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <algorithm>
+#include <map>
+#include <cassert>
/***********************************************************************
/
/ MACRO DEFINITIONS AND PARAMETERS
@@ -517,7 +524,6 @@
#define MBH_THERMAL 7
#define MBH_JETS 8
#define COLOR_FIELD 9
-#define SUPERNOVA_SEEDFIELD 11

/* Sink particle accretion modes */


diff -r bc76c0b9ebf0 -r 4c537ac40bce src/enzo/phys_constants.h
--- a/src/enzo/phys_constants.h
+++ b/src/enzo/phys_constants.h
@@ -60,4 +60,7 @@
#define kpc 3.0857e21
#define pc 3.0857e18

+/* Time units */
+#define SecondsPerYear 3.1556952e7
+
#endif

Repository URL: https://bitbucket.org/enzo/enzo-dev/

--

This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.