Commit 436c80b6 by Tom Theuns

parent 888a3fe5
 ... ... @@ -8,7 +8,8 @@ iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare f @physunits indir = './' basefile = 'Disk-Patch-dynamic_' ;basefile = 'Disc-Patch-dynamic_' basefile = 'Disc-Patch_' ; set properties of potential uL = phys.pc ; unit of length ... ... @@ -16,18 +17,27 @@ uM = phys.msun ; unit of mass uV = 1d5 ; unit of velocity ; properties of patch surface_density = 10. surface_density = 100. ; surface density of all mass, which generates the gravitational potential scale_height = 100. z_disk = 200.; z_disk = 200. ; fgas = 0.1 ; gas fraction gamma = 5./3. ; derived units constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ; pcentre = [0.,0.,z_disk] * pc / uL utherm = !pi * constG * surface_density * scale_height / (gamma-1.) temp = (utherm*uV^2)*phys.m_h/phys.kb soundspeed = sqrt(gamma * (gamma-1.) * utherm) t_dyn = sqrt(scale_height / (constG * surface_density)) rho0 = fgas*(surface_density)/(2.*scale_height) print,' dynamical time = ',t_dyn,' = ',t_dyn*UL/uV/(1d6*phys.yr),' Myr' print,' thermal energy per unit mass = ',utherm print,' central density = ',rho0,' = ',rho0*uM/uL^3/m_h,' particles/cm^3' print,' central temperature = ',temp lambda = 2 * !pi * phys.G^1.5 * (scale_height*uL)^1.5 * (surface_density * uM/uL^2)^0.5 * phys.m_h^2 / (gamma-1) / fgas print,' lambda = ',lambda stop ; infile = indir + basefile + '*' spawn,'ls -1 '+infile,res ... ...
 ... ... @@ -56,9 +56,10 @@ print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs # parameters of potential surface_density = 10. surface_density = 100. # surface density of all mass, which generates the gravitational potential scale_height = 100. gamma = 5./3. fgas = 0.1 # gas fraction # derived units const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs) ... ... @@ -131,7 +132,7 @@ h = glass_h[0:numGas] numGas = numpy.shape(pos)[0] # compute furthe properties of ICs column_density = surface_density * numpy.tanh(boxSize/2./scale_height) column_density = fgas * surface_density * numpy.tanh(boxSize/2./scale_height) enclosed_mass = column_density * boxSize * boxSize pmass = enclosed_mass / numGas meanrho = enclosed_mass / boxSize**3 ... ...
 ... ... @@ -97,7 +97,8 @@ //#define EXTERNAL_POTENTIAL_DISC_PATCH /* Source terms */ #define SN_FEEDBACK #define SOURCETERMS_NONE //#define SOURCETERMS_SN_FEEDBACK /* Cooling properties */ #define COOLING_NONE ... ...
 ... ... @@ -44,7 +44,7 @@ #include "potential.h" #include "runner.h" #include "scheduler.h" #include "sourceterms.h" #include "sourceterms_struct.h" #include "space.h" #include "task.h" #include "units.h" ... ...
 ... ... @@ -120,19 +120,11 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, * @param timer 1 if the time is to be recorded. */ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) { struct part *restrict parts = c->parts; struct xpart *restrict xparts = c->xparts; const int count = c->count; const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]}; const double cell_width[3] = {c->width[0], c->width[1], c->width[2]}; const int ti_current = r->e->ti_current; struct sourceterms *sourceterms = r->e->sourceterms; const double location[3] = {sourceterms->supernova.x, sourceterms->supernova.y, sourceterms->supernova.z}; const int dimen = 3; const double timeBase = r->e->timeBase; TIMER_TIC; ... ... @@ -143,78 +135,13 @@ void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) { return; } /* is supernova still active? */ if (sourceterms->supernova.status == supernova_is_not_done) { /* does cell contain explosion? */ if (count > 0) { const int incell = is_in_cell(cell_min, cell_width, location, dimen); if (incell == 1) { /* inject SN energy into particle with highest id, if it is active */ int imax = 0; struct part *restrict p_sn = NULL; struct xpart *restrict xp_sn = NULL; for (int i = 0; i < count; i++) { /* Get a direct pointer on the part. */ struct part *restrict p = &parts[i]; if (p->id > imax) { imax = p->id; p_sn = p; xp_sn = &xparts[i]; } } if (count > 0) { /* Is this part within the time step? */ if (p_sn->ti_begin == ti_current) { /* Does this time step straddle the feedback injection time? */ const float t_begin = p_sn->ti_begin * timeBase; const float t_end = p_sn->ti_end * timeBase; if (t_begin <= sourceterms->supernova.time && t_end > sourceterms->supernova.time) { /* store old time step */ const int dti_old = p_sn->ti_end - p_sn->ti_begin; /* add supernova feedback */ do_supernova_feedback(sourceterms, p_sn); /* label supernova as done */ sourceterms->supernova.status = supernova_is_done; message(" applied super nova, time = %d, location= %e %e %e", ti_current, p_sn->x[0], p_sn->x[1], p_sn->x[2]); message(" applied super nova, velocity = %e %e %e", p_sn->v[0], p_sn->v[1], p_sn->v[2]); /* update timestep if new time step shorter than old time step */ const int dti = get_part_timestep(p_sn, xp_sn, r->e); if (dti < dti_old) { p_sn->ti_end = p_sn->ti_begin + dti; message(" changed timestep from %d to %d", dti_old, dti); /* apply simple time-step limiter on all particles in same cell: */ int i_limit = 0; for (int i = 0; i < count; i++) { struct part *restrict p = &parts[i]; const int dti_old = p->ti_end - p->ti_begin; if (dti_old > 2 * dti) { i_limit++; const int dti_new = 2 * dti; p->ti_end = p->ti_begin + dti_new; message(" old step = %d new step = %d", dti_old, dti_new); } else message(" old step = %d", dti_old); } message(" count= %d limited timestep of %d particles ", count, i_limit); } // error(" check! "); } } } /* do sourceterms in this cell? */ const int incell = sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen); if (incell == 1) { sourceterms_apply(r, sourceterms, c); } } ... ...
 ... ... @@ -22,6 +22,7 @@ /* Local includes. */ #include "const.h" #include "hydro.h" #include "parser.h" #include "units.h" ... ... @@ -29,7 +30,7 @@ #include "sourceterms.h" /** * @brief Initialises the source terms * @brief Initialises the sourceterms * * @param parameter_file The parsed parameter file * @param us The current internal system of units ... ... @@ -37,30 +38,23 @@ */ void sourceterms_init(const struct swift_params* parameter_file, struct UnitSystem* us, struct sourceterms* source) { #ifdef SN_FEEDBACK source->supernova.time = parser_get_param_double(parameter_file, "SN:time"); source->supernova.energy = parser_get_param_double(parameter_file, "SN:energy"); source->supernova.x = parser_get_param_double(parameter_file, "SN:x"); source->supernova.y = parser_get_param_double(parameter_file, "SN:y"); source->supernova.z = parser_get_param_double(parameter_file, "SN:z"); source->supernova.status = supernova_is_not_done; #endif /* SN_FEEDBCK */ #ifdef SOURCETERMS_SN_FEEDBACK supernova_init(parameter_file, us, source); #endif /* SOURCETERMS_SN_FEEDBACK */ }; /** * @brief Prints the properties of the external potential to stdout. * * @brief Prints the properties of the source terms to stdout * @param source the structure that has all the source term properties */ void sourceterms_print(struct sourceterms* source) { #ifdef SN_FEEDBACK message( " Single SNe of energy= %e will explode at time= %e at location " "(%e,%e,%e)", source->supernova.energy, source->supernova.time, source->supernova.x, source->supernova.y, source->supernova.z); #endif /* SN_FEEDBACK */ #ifdef SOURCETERMS_NONE error(" no sourceterms defined yet you ran with -F"); #ifdef SOURCETERMS_SN_FEEDBACK #error can't have sourceterms when defined SOURCETERMS_NONE #endif #endif #ifdef SOURCETERMS_SN_FEEDBACK supernova_print(source); #endif /* SOURCETERMS_SN_FEEDBACK */ };
 ... ... @@ -18,22 +18,57 @@ ******************************************************************************/ #ifndef SWIFT_SOURCETERMS_H #define SWIFT_SOURCETERMS_H /** * @file src/sourceterms.h * @brief Branches between the different sourceterms functions. */ #include "./const.h" #ifdef SN_FEEDBACK #include "runner.h" #ifdef SOURCETERMS_SN_FEEDBACK #include "sourceterms/sn_feedback/sn_feedback_struct.h" #endif /* So far only one model here */ struct sourceterms { #ifdef SN_FEEDBACK #ifdef SOURCETERMS_SN_FEEDBACK struct supernova_struct supernova; #endif }; #ifdef SOURCETERMS_SN_FEEDBACK #include "sourceterms/sn_feedback/sn_feedback.h" #endif void sourceterms_init(const struct swift_params* parameter_file, struct UnitSystem* us, struct sourceterms* source); void sourceterms_print(struct sourceterms* source); #ifdef SN_FEEDBACK #include "sourceterms/sn_feedback/sn_feedback.h" /** * @file src/sourceterm.h * @brief Routines related to source terms * @param cell_min: corner of cell to test * @param cell_width: width of cell to test * @param sourceterms: properties of source terms to test * @param dimen: dimensionality of the problem * * This routine tests whether a source term should be applied to this cell * return: 1 if yes, return: 0 if no */ __attribute__((always_inline)) INLINE static int sourceterms_test_cell( const double cell_min[], const double cell_width[], struct sourceterms* sourceterms, const int dimen) { #ifdef SOURCETERMS_SN_FEEDBACK return supernova_feedback_test_cell(cell_min, cell_width, sourceterms, dimen); #endif return 0; }; __attribute__((always_inline)) INLINE static void sourceterms_apply( struct runner* r, struct sourceterms* sourceterms, struct cell* c) { #ifdef SOURCETERMS_SN_FEEDBACK supernova_feedback_apply(r, sourceterms, c); #endif }; #endif /* SWIFT_SOURCETERMS_H */
 /******************************************************************************* * This file is part of SWIFT. * Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk) * 2016 Tom Theuns (tom.theuns@durham.ac.uk) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with this program. If not, see . * ******************************************************************************/ #include #include "feedback.h" /** * @brief Computes the feedback time-step of a given particle due to * a source term * * This function only branches towards the potential chosen by the user. * * @param feedback The properties of the source terms. * @param phys_const The physical constants in internal units. * @param g Pointer to the particle data. */ __attribute__((always_inline)) INLINE static float sourceterms_compute_timestep( const struct sourceterms* feedback, const struct phys_const* const phys_const, const struct part* const p) { float dt = FLT_MAX; #ifdef SN_FEEDBACK dt = fmin(dt, sn_feedback_timestep(feedback, phys_const, p)); #endif } __attribute__((always_inline)) INLINE static void feedback( const struct sourceterms* feedback, const struct phys_const* const phys_const, struct part* p) #ifdef SN_FEEDBACK sn_feedback(feedback, phys_const, p); #endif }
 ... ... @@ -19,40 +19,174 @@ #ifndef SWIFT_SN_FEEDBACK_H #define SWIFT_SN_FEEDBACK_H #include /* Config parameters. */ #include "../config.h" #include "engine.h" #include "equation_of_state.h" #include "hydro.h" /* determine whether location is in cell */ __attribute__((always_inline)) INLINE static int is_in_cell( const double cell_min[], const double cell_width[], const double location[], const int dimen) { #include "runner.h" #include "timestep.h" /** * @file src/sourceterms/sn_feedback.h * * @brief Routines related to sourceterms (supernova feedback): determine if * feedback occurs in this cell * * @param cell_min: corner of cell to test * @param cell_width: width of cell to test * @param sourceterms: properties of source terms to test * @param dimen: dimensionality of the problem * * This routine tests whether a source term should be applied to this cell * return: 1 if yes, return: 0 if no */ __attribute__((always_inline)) INLINE static int supernova_feedback_test_cell( const double cell_min[], const double cell_width[], struct sourceterms* sourceterms, const int dimen) { if (sourceterms->supernova.status == supernova_is_done) return 0; const double location[3] = {sourceterms->supernova.x, sourceterms->supernova.y, sourceterms->supernova.z}; for (int i = 0; i < dimen; i++) { if (cell_min[i] > location[i]) return 0; if ((cell_min[i] + cell_width[i]) <= location[i]) return 0; } }; return 1; }; /** * @file src/sourceterms/sn_feedback.h * @brief Routines related to source terms (supernova feedback) * * @brief Routines related to source terms (supernova feedback): perform * feedback in this cell * @param r: the runner * @param sourceterms the structure describing the source terms properties * @param p the particle to apply feedback to * @param c the cell to apply feedback to * * This routine heats an individual particle (p), increasing its thermal energy * per unit mass * by supernova energy / particle mass. */ __attribute__((always_inline)) INLINE static void do_supernova_feedback( const struct sourceterms* sourceterms, struct part* p) { const float u_old = hydro_get_internal_energy(p, 0); message(" u_old= %e entropy= %e", u_old, p->entropy); const float u_new = u_old + sourceterms->supernova.energy / hydro_get_mass(p); hydro_set_internal_energy(p, u_new); const float u_set = hydro_get_internal_energy(p, 0.0); message(" unew = %e %e s= %e", u_new, u_set, p->entropy); message( " injected SN energy in particle = %lld, increased energy from %e to %e, " "check= %e", p->id, u_old, u_new, u_set); __attribute__((always_inline)) INLINE static void supernova_feedback_apply( struct runner* restrict r, struct sourceterms* restrict sourceterms, struct cell* restrict c) { const int count = c->count; struct part* restrict parts = c->parts; struct xpart* restrict xparts = c->xparts; const double timeBase = r->e->timeBase; const int ti_current = r->e->ti_current; /* inject SN energy into the particle with highest id in this cell if it is * active */ int imax = 0; struct part* restrict p_sn = NULL; struct xpart* restrict xp_sn = NULL; for (int i = 0; i < count; i++) { /* Get a direct pointer on the part. */ struct part* restrict p = &parts[i]; if (p->id > imax) { imax = p->id; p_sn = p; xp_sn = &xparts[i]; } } /* Is this part within the time step? */ if (p_sn->ti_begin == ti_current) { /* Does this time step straddle the feedback injection time? */ const float t_begin = p_sn->ti_begin * timeBase; const float t_end = p_sn->ti_end * timeBase; if (t_begin <= sourceterms->supernova.time && t_end > sourceterms->supernova.time) { /* store old time step */ const int dti_old = p_sn->ti_end - p_sn->ti_begin; /* add supernova feedback */ const float u_old = hydro_get_internal_energy(p_sn, 0); message(" u_old= %e entropy= %e", u_old, p_sn->entropy); const float u_new = u_old + sourceterms->supernova.energy / hydro_get_mass(p_sn); hydro_set_internal_energy(p_sn, u_new); const float u_set = hydro_get_internal_energy(p_sn, 0.0); const float ent_set = hydro_get_entropy(p_sn, 0.0); message(" unew = %e %e s= %e", u_new, u_set, ent_set); message( " injected SN energy in particle = %lld, increased energy from %e to " "%e, " "check= %e", p_sn->id, u_old, u_new, u_set); /* label supernova as done */ sourceterms->supernova.status = supernova_is_done; message(" applied super nova, time = %d, location= %e %e %e", ti_current, p_sn->x[0], p_sn->x[1], p_sn->x[2]); message(" applied super nova, velocity = %e %e %e", p_sn->v[0], p_sn->v[1], p_sn->v[2]); /* update timestep if new time step shorter than old time step */ const int dti = get_part_timestep(p_sn, xp_sn, r->e); if (dti < dti_old) { p_sn->ti_end = p_sn->ti_begin + dti; message(" changed timestep from %d to %d", dti_old, dti); /* apply simple time-step limiter on all particles in same cell: */ int i_limit = 0; for (int i = 0; i < count; i++) { struct part* restrict p = &parts[i]; const int dti_old = p->ti_end - p->ti_begin; if (dti_old > 2 * dti) { i_limit++; const int dti_new = 2 * dti; p->ti_end = p->ti_begin + dti_new; message(" old step = %d new step = %d", dti_old, dti_new); } else message(" old step = %d", dti_old); } message(" count= %d limited timestep of %d particles ", count, i_limit); } /* end of limiter */ } } }; /** * @file src/sourceterms/sn_feedback.h * * @brief Routine to initialise supernova feedback * @param parameterfile: the parse parmeter file * @param us: the unit system in use * @param sourceterms the structure describing the source terms properties * * This routine heats an individual particle (p), increasing its thermal energy * per unit mass * by supernova energy / particle mass. */ __attribute__((always_inline)) INLINE static void supernova_init( const struct swift_params* parameter_file, struct UnitSystem* us, struct sourceterms* source) { source->supernova.time = parser_get_param_double(parameter_file, "SN:time"); source->supernova.energy = parser_get_param_double(parameter_file, "SN:energy"); source->supernova.x = parser_get_param_double(parameter_file, "SN:x"); source->supernova.y = parser_get_param_double(parameter_file, "SN:y"); source->supernova.z = parser_get_param_double(parameter_file, "SN:z"); source->supernova.status = supernova_is_not_done; } __attribute__((always_inline)) INLINE static void supernova_print( struct sourceterms* source) { message( " Single SNe of energy= %e will explode at time= %e at location " "(%e,%e,%e)", source->supernova.energy, source->supernova.time, source->supernova.x, source->supernova.y, source->supernova.z); } #endif /* SWIFT_SN_FEEDBACK_H */
 ... ... @@ -24,8 +24,8 @@ /* Local headers. */ #include "const.h" #include "cooling.h" #include "debug.h" /** * @brief Compute a valid integer time-step form a given time-step * ... ...
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!