Commit 5380a052 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added a set of routines to manage unit conversion when doing i/o. Not used yet anywhere.


Former-commit-id: 3e7ffe02e2b791b3f07c014515d0fff43ad751f5
parent 2abe8136
......@@ -714,6 +714,7 @@ int main ( int argc , char *argv[] ) {
struct part *parts = NULL;
struct space s;
struct engine e;
struct UnitSystem us;
char ICfileName[200];
float dt_max = 0.0f;
ticks tic;
......@@ -826,6 +827,18 @@ int main ( int argc , char *argv[] ) {
if ( myrank == 0 )
message( "sizeof(struct part) is %li bytes." , (long int)sizeof( struct part ));
/* Initilaize unit system */
initUnitSystem(&us);
if ( myrank == 0 )
{
message( "Unit system: U_L = %e cm.", us.UnitLength_in_cgs );
message( "Unit system: U_M = %e g.", us.UnitMass_in_cgs );
message( "Unit system: U_t = %e s.", us.UnitTime_in_cgs );
message( "Unit system: U_I = %e A.", us.UnitCurrent_in_cgs );
message( "Unit system: U_T = %e K.", us.UnitTemperature_in_cgs );
/* message( "Density units: %e a^%f h^%f.", conversionFactor(&us, UNIT_CONV_ENTROPY), aFactor(&us, UNIT_CONV_ENTROPY), hFactor(&us, UNIT_CONV_ENTROPY) ); */
}
/* Read particles and space information from (GADGET) IC */
tic = getticks();
#ifdef WITH_MPI
......@@ -840,7 +853,7 @@ int main ( int argc , char *argv[] ) {
if( scaling != 1.0 )
for ( k = 0 ; k < N ; k++ )
parts[k].h *= scaling;
/* Apply shift */
if(shift[0] !=0 || shift[1] !=0 || shift[2] !=0 )
for ( k = 0 ; k < N ; k++ ) {
......
......@@ -42,7 +42,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
io.c timers.c debug.c scheduler.c proxy.c parallel_io.c
io.c timers.c debug.c scheduler.c proxy.c parallel_io.c units.c
# Sources and flags for regular library
libswiftsim_la_SOURCES = $(AM_SOURCES)
......
......@@ -43,3 +43,9 @@
/* SPH variant to use */
#define LEGACY_GADGET2_SPH
/* System of units */
#define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */
#define const_unit_mass_in_cgs 1 /* 1.9891e33 /\* 1 M_sun *\/ */
#define const_unit_velocity_in_cgs 1 /* 1e5 /\* km s^-1 *\/ */
......@@ -39,6 +39,7 @@
#include "io.h"
#include "parallel_io.h"
#include "debug.h"
#include "units.h"
#ifdef LEGACY_GADGET2_SPH
#include "runner_iact_legacy.h"
......
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
* Matthieu Schaller (matthieu.schaller@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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
/* Config parameters. */
#include "../config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stddef.h>
#include <math.h>
#include "const.h"
#include "cycle.h"
#include "part.h"
#include "error.h"
#include "units.h"
/**
* @brief Initialises the UnitSystem structure with the constants given in const.h
* @param us The UnitSystem to initialize
*/
void initUnitSystem(struct UnitSystem* us)
{
us->UnitMass_in_cgs = const_unit_mass_in_cgs;
us->UnitLength_in_cgs = const_unit_length_in_cgs;
us->UnitTime_in_cgs = 1.d / (const_unit_velocity_in_cgs / const_unit_length_in_cgs );
us->UnitCurrent_in_cgs = 1.;
us->UnitTemperature_in_cgs = 1.;
}
void getBaseUnitExponantsArray(float baseUnitsExp[5], enum UnitConversionFactor unit)
{
switch( unit )
{
case UNIT_CONV_MASS:
baseUnitsExp[UNIT_MASS] = 1; break;
case UNIT_CONV_LENGTH:
baseUnitsExp[UNIT_LENGTH] = 1; break;
case UNIT_CONV_TIME:
baseUnitsExp[UNIT_TIME] = 1; break;
case UNIT_CONV_FREQUENCY:
baseUnitsExp[UNIT_TIME] = -1; break;
case UNIT_CONV_DENSITY:
baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = -3; break;
case UNIT_CONV_SPEED:
baseUnitsExp[UNIT_LENGTH] = 1; baseUnitsExp[UNIT_TIME] = -1; break;
case UNIT_CONV_ACCELERATION:
baseUnitsExp[UNIT_LENGTH] = 1; baseUnitsExp[UNIT_TIME] = -2; break;
case UNIT_CONV_FORCE:
baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = 1; baseUnitsExp[UNIT_TIME] = -2; break;
case UNIT_CONV_ENERGY:
baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = 2; baseUnitsExp[UNIT_TIME] = -2; break;
case UNIT_CONV_ENERGY_PER_UNIT_MASS:
baseUnitsExp[UNIT_LENGTH] = 2; baseUnitsExp[UNIT_TIME] = -2; break;
case UNIT_CONV_ENTROPY:
baseUnitsExp[UNIT_MASS] = 1. - const_hydro_gamma; baseUnitsExp[UNIT_LENGTH] = 3.*const_hydro_gamma - 1.; baseUnitsExp[UNIT_TIME] = -2; break;
case UNIT_CONV_POWER:
baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = 2; baseUnitsExp[UNIT_TIME] = -3; break;
case UNIT_CONV_PRESSURE:
baseUnitsExp[UNIT_MASS] = 1; baseUnitsExp[UNIT_LENGTH] = -1; baseUnitsExp[UNIT_TIME] = -2; break;
}
}
/**
* @brief Returns the conversion factor for a given unit in the chosen unit system
* @param us The system of units in use
* @param unit The unit to convert
*/
double conversionFactor(struct UnitSystem* us, enum UnitConversionFactor unit)
{
float baseUnitsExp[5] = { 0, 0, 0, 0, 0 };
getBaseUnitExponantsArray(baseUnitsExp, unit);
return generalConversionFactor(us, baseUnitsExp);
}
/**
* @brief Returns the h factor exponentiation for a given unit
* @param us The system of units in use
* @param unit The unit to convert
*/
double hFactor(struct UnitSystem* us, enum UnitConversionFactor unit)
{
float baseUnitsExp[5] = { 0, 0, 0, 0, 0 };
getBaseUnitExponantsArray(baseUnitsExp, unit);
return generalhFactor(us, baseUnitsExp);
}
/**
* @brief Returns the scaling factor exponentiation for a given unit
* @param us The system of units in use
* @param unit The unit to convert
*/
double aFactor(struct UnitSystem* us, enum UnitConversionFactor unit)
{
float baseUnitsExp[5] = { 0, 0, 0, 0, 0 };
getBaseUnitExponantsArray(baseUnitsExp, unit);
return generalaFactor(us, baseUnitsExp);
}
/**
* @brief Returns the conversion factor for a given unit (expressed in terms of the 5 fundamental units) in the chosen unit system
* @param us The unit system used
* @param baseUnitsExponants The exponant (integer) of each base units required to form the desired quantity. See conversionFactor() for a working example
*/
double generalConversionFactor(struct UnitSystem* us, float baseUnitsExponants[5])
{
double factor = 1.;
int i;
for(i = 0 ; i < 5 ; ++i )
if(baseUnitsExponants[i] != 0)
factor *= pow( *( &( us->UnitMass_in_cgs ) + i ) , baseUnitsExponants[i] );
return factor;
}
/**
* @brief Returns the h factor exponentiation for a given unit (expressed in terms of the 5 fundamental units)
* @param us The unit system used
* @param baseUnitsExponants The exponant (integer) of each base units required to form the desired quantity. See conversionFactor() for a working example
*/
double generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5])
{
float factor_exp = 0;
factor_exp += -baseUnitsExponants[UNIT_MASS];
factor_exp += -baseUnitsExponants[UNIT_LENGTH];
factor_exp += -baseUnitsExponants[UNIT_TIME];
return factor_exp;
}
/**
* @brief Returns the scaling factor exponentiation for a given unit (expressed in terms of the 5 fundamental units)
* @param us The unit system used
* @param baseUnitsExponants The exponant (integer) of each base units required to form the desired quantity. See conversionFactor() for a working example
*/
double generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5])
{
float factor_exp = 0;
factor_exp += baseUnitsExponants[UNIT_LENGTH];
return factor_exp;
}
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2012 Matthieu Schaller (matthieu.schaller@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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
/**
* @brief The unit system used internally.
*
* This structure contains the conversion factors to the 7 cgs base units to the internal units.
* It is used everytime a conversion is performed or an i/o function is called.
*
**/
struct UnitSystem
{
double UnitMass_in_cgs; /*< Conversion factor from grams to internal mass units */
double UnitLength_in_cgs; /*< Conversion factor from centimeters to internal length units. */
double UnitTime_in_cgs; /*< Conversion factor from seconds to internal time units. */
double UnitCurrent_in_cgs; /*< Conversion factor from Ampere to internal current units. */
double UnitTemperature_in_cgs; /*< Conversion factor from Kelvins to internal temperature units. */
};
/**
* @brief The base units used in the cgs (and internal) system. All units are derived from those.
*/
enum BaseUnits
{
UNIT_MASS,
UNIT_LENGTH,
UNIT_TIME,
UNIT_CURRENT,
UNIT_TEMPERATURE
};
/**
* @brief The different conversion factors supported by default
*/
enum UnitConversionFactor
{
UNIT_CONV_MASS,
UNIT_CONV_LENGTH,
UNIT_CONV_TIME,
UNIT_CONV_DENSITY,
UNIT_CONV_SPEED,
UNIT_CONV_ACCELERATION,
UNIT_CONV_FORCE,
UNIT_CONV_ENERGY,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
UNIT_CONV_ENTROPY,
UNIT_CONV_POWER,
UNIT_CONV_PRESSURE,
UNIT_CONV_FREQUENCY
};
/**
* @brief Initialises the UnitSystem structure with the constants given in const.h
*/
void initUnitSystem(struct UnitSystem*);
/**
* @brief Returns the conversion factor for a given unit (expressed in terms of the 5 fundamental units) in the chosen unit system
*/
double generalConversionFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
/**
* @brief Returns the conversion factor for a given unit in the chosen unit system
*/
double conversionFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
/**
* @brief Returns the h factor for a given unit (expressed in terms of the 5 fundamental units) in the chosen unit system
*/
double generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
/**
* @brief Returns the h factor for a given unit in the chosen unit system
*/
double hFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
/**
* @brief Returns the scaling factor for a given unit (expressed in terms of the 5 fundamental units) in the chosen unit system
*/
double generalaFactor(struct UnitSystem* us, float baseUnitsExponants[5]);
/**
* @brief Returns the scaling factor for a given unit in the chosen unit system
*/
double aFactor(struct UnitSystem* us, enum UnitConversionFactor unit);
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment