[cig-commits] [commit] master: 3D support (12cc2dd)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue May 20 12:19:53 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/a21aaa79d75b89a3d50d5c865e7dc51d309da9b9...b2eaff9e459f8351633e8b4b43c1284c90373873
>---------------------------------------------------------------
commit 12cc2dd8aeb0093328940c69876995f3b6f53c83
Author: Siqi Zhang <siqi.zhang at mq.edu.au>
Date: Thu May 15 17:40:48 2014 -0500
3D support
>---------------------------------------------------------------
12cc2dd8aeb0093328940c69876995f3b6f53c83
source/initial_conditions/solidus.cc | 198 ++++++++++++++++++-----------------
1 file changed, 101 insertions(+), 97 deletions(-)
diff --git a/source/initial_conditions/solidus.cc b/source/initial_conditions/solidus.cc
index e69d913..8dc2c3f 100644
--- a/source/initial_conditions/solidus.cc
+++ b/source/initial_conditions/solidus.cc
@@ -27,64 +27,64 @@
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);
- char temp[256];
- std::string T_Unit,P_Unit;
+ 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);
+ char temp[256];
+ std::string T_Unit,P_Unit;
Num_points=0;
- if(in.fail())return;
- in.getline(temp,256);
- in>>T_Unit>>P_Unit;
- in.getline(temp,256);
+ if(in.fail())return;
+ in.getline(temp,256);
+ in>>T_Unit>>P_Unit;
+ in.getline(temp,256);
while(!in.eof())
{
- double T,p;
+ 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")
- {
- is_radius=true; // Second column in radius instead of pressure
- p*=1.e3; // km to meters
- }
- T_array.push_back(T);
- P_array.push_back(p);
- Num_points++;
+ //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")
+ {
+ is_radius=true; // Second column in radius instead of pressure
+ p*=1.e3; // km to meters
+ }
+ T_array.push_back(T);
+ P_array.push_back(p);
+ Num_points++;
}
- in.getline(temp,256);
+ 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.);
- for(unsigned i=1;i<Num_points;i++)
- {
- if( (i==Num_points-1) ||
- (is_radius && P_value>P_array[i]) ||
- (!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);
- }
-
- }
+ 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.);
+ for(unsigned i=1;i<Num_points;i++)
+ {
+ if( (i==Num_points-1) ||
+ (is_radius && P_value>P_array[i]) ||
+ (!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);
+ }
+
+ }
namespace InitialConditions
{
@@ -97,30 +97,34 @@ namespace aspect
Solidus<dim>::
initial_temperature (const Point<dim> &position) const
{
- double R1,T_min,T_litho;
- double Theta=atan2(position(0),position(1));
- double T_solidus,T_perturbation;
- double litho_thick_theta;
- 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."));
+ double R1,T_min,T_litho;
+ double Theta=atan2(position(0),position(1));
+ double Phi;
+ if(dim==3)
+ Phi=atan2(position(0),position(2));
+ else
+ Phi=atan2(1,0);
+ double T_solidus,T_perturbation;
+ double litho_thick_theta;
+ 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(dynamic_cast< const GeometryModel::SphericalShell<dim> *>( &this->get_geometry_model() )!=0,
- ExcMessage("This initial condition can only be work with sphereical shell geometry model."));
+ ExcMessage("This initial condition can only be work with sphereical shell geometry model."));
R1=(dynamic_cast< const GeometryModel::SphericalShell<dim> &>(this->get_geometry_model())).R1;
- T_min=(this->get_boundary_temperature()).minimal_temperature();
+ T_min=(this->get_boundary_temperature()).minimal_temperature();
- litho_thick_theta=litho_thick-Magnitude_lith*sin(n*Theta);
- T_litho=Solidus_curve.T(0,R1-litho_thick_theta)+deltaT;
+ litho_thick_theta=litho_thick-Magnitude_lith*sin(n*Theta);
+ T_litho=Solidus_curve.T(0,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;
+ 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;
- //Currently only work in 2D
- T_perturbation=Depth/( this->geometry_model->maximal_depth() )*Magnitude_T*sin(n*Theta);
+ T_perturbation=Depth/( this->geometry_model->maximal_depth() )*Magnitude_T*sin(n*Theta)*sin(n*Phi);
return T_solidus+T_perturbation;
}
@@ -130,35 +134,35 @@ namespace aspect
{
prm.enter_subsection("Initial conditions");
{
- prm.declare_entry ("Supersolidus","0e0",
- Patterns::Double (),
- "The difference from solidus.");
+ prm.declare_entry ("Supersolidus","0e0",
+ Patterns::Double (),
+ "The difference from solidus.");
prm.declare_entry ("Lithosphere thickness","0",
Patterns::Double (0),
"The thickness of lithosphere thickness. Unit: m");
- prm.enter_subsection("Perturbation");
+ prm.enter_subsection("Perturbation");
{
- prm.declare_entry ("Temperature amplitude", "0e0",
- Patterns::Double (0),
- "The amplitude of the initial spherical temperature perturbation in (K)");
- prm.declare_entry ("Lithosphere thickness amplitude", "0e0",
- Patterns::Double (),
- "The amplitude of the initial lithosphere thickness perturbation in (m)");
- prm.declare_entry ("Degree","1",
- Patterns::Integer(0),
- "The degree of the perturbation.");
+ prm.declare_entry ("Temperature amplitude", "0e0",
+ Patterns::Double (0),
+ "The amplitude of the initial spherical temperature perturbation in (K)");
+ prm.declare_entry ("Lithosphere thickness amplitude", "0e0",
+ Patterns::Double (),
+ "The amplitude of the initial lithosphere thickness perturbation in (m)");
+ prm.declare_entry ("Degree","1",
+ Patterns::Integer(0),
+ "The degree of the perturbation.");
}
prm.leave_subsection();
- prm.enter_subsection ("Data");
- {
- prm.declare_entry ("Solidus filename", "",
+ prm.enter_subsection ("Data");
+ {
+ prm.declare_entry ("Solidus filename", "",
Patterns::Anything(),
"The solidus filename.");
}
prm.leave_subsection();
}
- prm.leave_subsection();
- }
+ prm.leave_subsection();
+ }
template <int dim>
@@ -167,19 +171,19 @@ namespace aspect
{
prm.enter_subsection("Initial conditions");
{
- deltaT=prm.get_double("Supersolidus");
- 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");
- n = prm.get_integer("Degree");
-
- }
- prm.leave_subsection();
+ deltaT=prm.get_double("Supersolidus");
+ 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");
+ n = prm.get_integer("Degree");
+
+ }
+ prm.leave_subsection();
prm.enter_subsection("Data");
{
- solidus_filename=prm.get ("Solidus filename");
+ solidus_filename=prm.get ("Solidus filename");
}
prm.leave_subsection();
}
@@ -197,6 +201,6 @@ namespace aspect
ASPECT_REGISTER_INITIAL_CONDITIONS(Solidus,
"solidus",
"Temperature initial condition as solidus,"
- "with perturbation as sin() funciton.");
+ "with perturbation as sin() funciton.");
}
}
More information about the CIG-COMMITS
mailing list