[cig-commits] [commit] master: bug fix (a9bf3ee)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue May 20 12:20:21 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/a21aaa79d75b89a3d50d5c865e7dc51d309da9b9...b2eaff9e459f8351633e8b4b43c1284c90373873
>---------------------------------------------------------------
commit a9bf3ee880ca471cb0f607d19c16ae38e37b3bb0
Author: Siqi Zhang <siqi.zhang at mq.edu.au>
Date: Tue May 20 10:59:23 2014 -0500
bug fix
>---------------------------------------------------------------
a9bf3ee880ca471cb0f607d19c16ae38e37b3bb0
include/aspect/initial_conditions/solidus.h | 116 +++++++++++++++++++++-------
source/initial_conditions/solidus.cc | 116 +++++++++++++++++-----------
2 files changed, 160 insertions(+), 72 deletions(-)
diff --git a/include/aspect/initial_conditions/solidus.h b/include/aspect/initial_conditions/solidus.h
index dd330fe..0d00c3b 100644
--- a/include/aspect/initial_conditions/solidus.h
+++ b/include/aspect/initial_conditions/solidus.h
@@ -28,12 +28,57 @@
namespace aspect
{
+
+
namespace InitialConditions
{
+ /**
+ * Data class to handle the melting curve.
+ */
+ class Melting_curve
+ {
+ public:
+
+ /**
+ * Read the data file into the class.
+ */
+ void read(const std::string &filename);
+
+ /**
+ * Get the melting temperature.
+ */
+ double T(const double p, const double radius) const;
+
+ /**
+ * Is the melting curve denpendent on radius.
+ * The melting curve can be dependent on radius or pressure.
+ */
+ bool is_radius;
+
+ /**
+ * Number of data points in the melting curve data.
+ */
+ unsigned int Num_points;
+ private:
+ /**
+ * Data array for temperature.
+ */
+ std::vector<double> T_array;
+
+ /**
+ * Data array for pressure/radius.
+ */
+ std::vector<double> P_array;
+
+ /**
+ * Name of the data file.
+ */
+ std::string data_filename;
+ };
+
/**
- * A class that implements temperature initial conditions based on a
- * functional description provided in the input file.
+ * A class that implements temperature initial conditions based on solidus provided by a data file.
*
* @ingroup InitialConditionsModels
*/
@@ -79,32 +124,51 @@ namespace aspect
const Tensor<1,dim>
spherical_surface_coordinates(const Tensor<1,dim> &position) const;
- double litho_thick;
- double Magnitude_T;
- double Magnitude_lith;
- double deltaT;
- int lateral_wave_number_1;
- int lateral_wave_number_2;
- std::string solidus_filename;
+ /**
+ * Lithosphere thickness.
+ */
+ double litho_thick;
+
+ /**
+ * Magnitude of temperature perturbation.
+ */
+ double magnitude_T;
+
+ /**
+ * Magnitude of lithosphere thickness perturbation.
+ */
+ double magnitude_lith;
+
+ /**
+ * Temperature difference from solidus, so the initial condition can be super-solidus or sub-solidus.
+ */
+ double deltaT;
+
+ /**
+ * The lateral wave number of the harmonic perturbation in the first
+ * dimension. This is the only lateral wave number in 2D and equals
+ * the degree of the spherical harmonics in a 3D spherical shell.
+ */
+ int lateral_wave_number_1;
+
+ /**
+ * The lateral wave number of the harmonic perturbation in the second
+ * dimension. This is not used in 2D and equals the order of the
+ * spherical harmonics in a 3D spherical shell.
+ */
+ int lateral_wave_number_2;
+
+ /**
+ * The data file name for solidus data.
+ */
+ std::string solidus_filename;
+
+ /**
+ * Data class for melting curve
+ */
+ Melting_curve solidus_curve;
};
}
- namespace melting
- {
- class Melting_curve
- {
- public:
- Melting_curve(const std::string &filename);
- void read(const std::string &filename);
- double T(const double p, const double radius) const;
- bool is_radius;
- unsigned int Num_points;
- private:
- //unsigned int Num_points;
- std::vector<double> T_array;
- std::vector<double> P_array;
- //bool is_radius;
- };
- }
}
diff --git a/source/initial_conditions/solidus.cc b/source/initial_conditions/solidus.cc
index b0b9ac5..117aaf9 100644
--- a/source/initial_conditions/solidus.cc
+++ b/source/initial_conditions/solidus.cc
@@ -30,48 +30,61 @@ namespace aspect
{
namespace melting
{
- Melting_curve::Melting_curve(const std::string &filename)
- {
- read(filename);
- }
void Melting_curve::read(const std::string &filename)
{
- std::ifstream in(filename.c_str(), std::ios::in);
+ data_filename=filename;
+ std::ifstream in(data_filename.c_str(), std::ios::in);
char temp[256];
std::string T_Unit,P_Unit;
- Num_points=0;
+ Num_points=0;
if(in.fail())return;
in.getline(temp,256);
in>>T_Unit>>P_Unit;
in.getline(temp,256);
- while(!in.eof())
- {
+ while(!in.eof())
+ {
double T,p;
in>>T>>p;
if(!in.fail())
{
//Unit switching
- if(T_Unit=="C") T+=273.15; // Degree C to K
- if(P_Unit=="kbar") p*=1.e8; // kbar to Pa
- if(P_Unit=="GPa") p*=1.e9; // GPa to Pa
- is_radius=false; // Second column is pressure
- if(P_Unit=="km")
+ if(T_Unit=="C")
{
- is_radius=true; // Second column in radius instead of pressure
- p*=1.e3; // km to meters
+ T+=273.15; // Degree C to K
}
+ else if(T_Unit!="K")
+ {
+ AssertThrow(false,ExcMessage ("Unit of the first column of melting curve data "
+ "has to be one of the following: C/K."))
+ }
+
+ is_radius=false; // Second column is pressure
+ if(P_Unit=="kbar")
+ p*=1.e8; // kbar to Pa
+ else if(P_Unit=="GPa")
+ p*=1.e9; // GPa to Pa
+ else
+ {
+ is_radius=true; // Second column in radius instead of pressure
+ if(P_Unit=="km")
+ p*=1.e3; // km to meters
+ else if(P_Unit!="m")
+ AssertThrow(false,ExcMessage ("Unit of the second column of melting curve data "
+ "has to be one of the following: Pa/Gpa/km/m."))
T_array.push_back(T);
P_array.push_back(p);
Num_points++;
- }
+ }
in.getline(temp,256);
- }
- }
+ }
+ }
+
+
double Melting_curve::T(const double p, const double radius) const
{
double T_value,P_value=is_radius?radius:p;
- if(T_array.size()==0)return(0.);
+ if(T_array.size()==0)return(0);
for(unsigned i=1;i<Num_points;i++)
{
if( (i==Num_points-1) ||
@@ -79,10 +92,11 @@ namespace aspect
(!is_radius && P_value<P_array[i]) )
{
T_value=T_array[i-1]+(T_array[i]-T_array[i-1])/(P_array[i]-P_array[i-1])*(P_value-P_array[i-1]);
- break;
+ return(T_value);
}
}
- return(T_value);
+ AssertThrow(false,ExcMessage(std::string("Something wrong with the melting curve data ")+ data_filename ));
+ return(-1);
}
}
@@ -91,6 +105,7 @@ namespace aspect
{
template <int dim>
Solidus<dim>::Solidus ()
+ solidus_curve()
{}
template <int dim>
@@ -101,32 +116,31 @@ namespace aspect
double T_min,T_litho;
double T_solidus,T_perturbation;
double litho_thick_theta;
- double lateral_perturbation;
+ double lateral_perturbation;
double Depth=this->geometry_model->depth(position);
- static melting::Melting_curve Solidus_curve(solidus_filename);
- AssertThrow(Solidus_curve.is_radius==true,ExcMessage("The solidus curve has to be depth dependent."));
- AssertThrow(Solidus_curve.Num_points!=0,ExcMessage("Error eading solidus file."));
+ AssertThrow(solidus_curve.is_radius==true,ExcMessage("The solidus curve has to be radius dependent."));
+ AssertThrow(solidus_curve.Num_points!=0,ExcMessage("Error eading solidus file."));
const GeometryModel::SphericalShell<dim> *spherical_geometry_model=
- dynamic_cast< const GeometryModel::SphericalShell<dim> *>(this->geometry_model);
- AssertThrow(spherical_geometry_model!=0,
+ dynamic_cast< const GeometryModel::SphericalShell<dim> *>(this->geometry_model);
+ AssertThrow(spherical_geometry_model!=0,
ExcMessage("This initial condition can only be work with spherical shell geometry model."));
T_min=(this->get_boundary_temperature()).minimal_temperature();
// In case of spherical shell calculate spherical coordinates
- const Tensor<1,dim> scoord = spherical_surface_coordinates(position);
- if(dim==2)
- {
+ const Tensor<1,dim> scoord = spherical_surface_coordinates(position);
+ if(dim==2)
+ {
// Use a sine as lateral perturbation that is scaled to the opening angle of the geometry.
// This way the perturbation is alway 0 at the model boundaries.
const double opening_angle = spherical_geometry_model->opening_angle()*numbers::PI/180.0;
lateral_perturbation = std::sin(lateral_wave_number_1*scoord[1]*numbers::PI/opening_angle);
- }
- else if(dim==3)
- {
- // Spherical harmonics are only defined for order <= degree
+ }
+ else if(dim==3)
+ {
+ // Spherical harmonics are only defined for order <= degree
// and degree >= 0. Verify that it is indeed.
- Assert ( std::abs(lateral_wave_number_2) <= lateral_wave_number_1,
+ Assert ( std::abs(lateral_wave_number_2) <= lateral_wave_number_1,
ExcMessage ("Spherical harmonics can only be computed for "
"order <= degree."));
Assert ( lateral_wave_number_1 >= 0,
@@ -134,18 +148,21 @@ namespace aspect
"degree >= 0."));
// use a spherical harmonic function as lateral perturbation
lateral_perturbation = boost::math::spherical_harmonic_r(lateral_wave_number_1,lateral_wave_number_2,scoord[2],scoord[1]);
- }
- litho_thick_theta=litho_thick-Magnitude_lith*lateral_perturbation;
- T_litho=Solidus_curve.T(0,spherical_geometry_model->R1-litho_thick_theta)+deltaT;
+ }
+ litho_thick_theta=litho_thick-magnitude_lith*lateral_perturbation;
+ T_litho=solidus_curve.T(0,spherical_geometry_model->R1-litho_thick_theta)+deltaT;
if(litho_thick_theta>0 && Depth<litho_thick_theta)
T_solidus=T_min+(T_litho-T_min)*(Depth/litho_thick_theta);
else
- T_solidus=Solidus_curve.T(0,sqrt(position.square()))+deltaT;
+ T_solidus=solidus_curve.T(0,sqrt(position.square()))+deltaT;
- T_perturbation=Depth/( this->geometry_model->maximal_depth() )*Magnitude_T*lateral_perturbation;
+ T_perturbation=Depth/( this->geometry_model->maximal_depth() )*magnitude_T*lateral_perturbation;
return T_solidus+T_perturbation;
}
+
+
+
template <int dim>
const Tensor<1,dim>
Solidus<dim>::
@@ -153,15 +170,17 @@ namespace aspect
{
Tensor<1,dim> scoord;
- scoord[0] = std::sqrt(position.norm_square()); // R
- scoord[1] = std::atan2(position[1],position[0]); // Phi
- if (scoord[1] < 0.0) scoord[1] = 2*numbers::PI + scoord[1]; // correct phi to [0,2*pi]
+ scoord[0] = std::sqrt(position.norm_square()); // R
+ scoord[1] = std::atan2(position[1],position[0]); // Phi
+ if (scoord[1] < 0.0) scoord[1] = 2*numbers::PI + scoord[1]; // correct phi to [0,2*pi]
if (dim==3)
scoord[2] = std::acos(position[2]/std::sqrt(position.norm_square())); // Theta
return scoord;
}
+
+
template <int dim>
void
Solidus<dim>::declare_parameters (ParameterHandler &prm)
@@ -204,13 +223,14 @@ namespace aspect
prm.declare_entry ("Solidus filename", "",
Patterns::Anything(),
"The solidus filename.");
- }
+ }
prm.leave_subsection();
}
prm.leave_subsection();
}
+
template <int dim>
void
Solidus<dim>::parse_parameters (ParameterHandler &prm)
@@ -221,8 +241,8 @@ namespace aspect
litho_thick=prm.get_double("Lithosphere thickness");
prm.enter_subsection("Perturbation");
{
- Magnitude_T = prm.get_double("Temperature amplitude");
- Magnitude_lith = prm.get_double("Lithosphere thickness amplitude");
+ magnitude_T = prm.get_double("Temperature amplitude");
+ magnitude_lith = prm.get_double("Lithosphere thickness amplitude");
lateral_wave_number_1 = prm.get_integer ("Lateral wave number one");
lateral_wave_number_2 = prm.get_integer ("Lateral wave number two");
}
@@ -230,6 +250,7 @@ namespace aspect
prm.enter_subsection("Data");
{
solidus_filename=prm.get ("Solidus filename");
+ solidus_curve.read(solidus_filename);
}
prm.leave_subsection();
}
@@ -239,6 +260,9 @@ namespace aspect
}
}
+
+
+
// explicit instantiations
namespace aspect
{
More information about the CIG-COMMITS
mailing list