[cig-commits] [commit] master: .. (f701f89)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue May 20 12:20:10 PDT 2014
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/a21aaa79d75b89a3d50d5c865e7dc51d309da9b9...b2eaff9e459f8351633e8b4b43c1284c90373873
>---------------------------------------------------------------
commit f701f89854137cea6c6653f9aca4076a4a47fac3
Author: Siqi Zhang <siqi.zhang at mq.edu.au>
Date: Mon May 19 12:25:08 2014 -0500
..
>---------------------------------------------------------------
f701f89854137cea6c6653f9aca4076a4a47fac3
source/initial_conditions/solidus.cc | 252 -----------------------------------
1 file changed, 252 deletions(-)
diff --git a/source/initial_conditions/solidus.cc b/source/initial_conditions/solidus.cc
deleted file mode 100644
index 50f039b..0000000
--- a/source/initial_conditions/solidus.cc
+++ /dev/null
@@ -1,252 +0,0 @@
-/*
- Copyright (C) 2011, 2012 by the authors of the ASPECT code.
-
- This file is part of ASPECT.
-
- ASPECT is free software; you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation; either version 2, or (at your option)
- any later version.
-
- ASPECT 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 General Public License
- along with ASPECT; see the file doc/COPYING. If not see
- <http://www.gnu.org/licenses/>.
-*/
-/* $Id: function.cc 1389 2012-11-28 14:08:16Z bangerth $ */
-
-//#define _USE_MATH_DEFINES
-#include <cmath>
-#include <aspect/initial_conditions/solidus.h>
-#include <aspect/geometry_model/spherical_shell.h>
-#include <aspect/boundary_temperature/interface.h>
-#include <boost/math/special_functions/spherical_harmonic.hpp>
-
-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;
- Num_points=0;
- if(in.fail())return;
- in.getline(temp,256);
- in>>T_Unit>>P_Unit;
- in.getline(temp,256);
- 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")
- {
- 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);
- }
- }
-
- 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
- {
- template <int dim>
- Solidus<dim>::Solidus ()
- {}
-
- template <int dim>
- double
- Solidus<dim>::
- initial_temperature (const Point<dim> &position) const
- {
- double T_min,T_litho;
- double T_solidus,T_perturbation;
- double litho_thick_theta;
- 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."));
- const GeometryModel::SphericalShell<dim> *spherical_geometry_model=
- dynamic_cast< const GeometryModel::SphericalShell<dim> *>(this->geometry_model);
- AssertThrow(spherical_geometry_model!=0,
- ExcMessage("This initial condition can only be work with sphereical 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)
- {
- // 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
- // and degree >= 0. Verify that it is indeed.
- 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,
- ExcMessage ("Spherical harmonics can only be computed for "
- "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;
-
- 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_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>::
- spherical_surface_coordinates(const Tensor<1,dim> &position) const
- {
- 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]
- 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)
- {
- prm.enter_subsection("Initial conditions");
- {
- 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.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 ("Lateral wave number one","3",
- Patterns::Integer(),
- "Doubled first lateral wave number of the harmonic perturbation. "
- "Equals the spherical harmonic degree in 3D spherical shells. "
- "In all other cases one equals half of a sine period over "
- "the model domain. This allows for single up-/downswings. "
- "Negative numbers reverse the sign of the perturbation but are "
- "not allowed for the spherical harmonic case.");
- prm.declare_entry ("Lateral wave number two", "2",
- Patterns::Integer (),
- "Doubled second lateral wave number of the harmonic perturbation. "
- "Equals the spherical harmonic order in 3D spherical shells. "
- "In all other cases one equals half of a sine period over "
- "the model domain. This allows for single up-/downswings. "
- "Negative numbers reverse the sign of the perturbation.");
- }
- prm.leave_subsection();
- prm.enter_subsection ("Data");
- {
- 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)
- {
- 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");
- lateral_wave_number_1 = prm.get_integer ("Lateral wave number one");
- lateral_wave_number_2 = prm.get_integer ("Lateral wave number two");
- }
- prm.leave_subsection();
- prm.enter_subsection("Data");
- {
- solidus_filename=prm.get ("Solidus filename");
- }
- prm.leave_subsection();
- }
- prm.leave_subsection();
- }
-
- }
-}
-
-// explicit instantiations
-namespace aspect
-{
- namespace InitialConditions
- {
- ASPECT_REGISTER_INITIAL_CONDITIONS(Solidus,
- "solidus",
- "Temperature initial condition as solidus,"
- "with perturbation as sin() funciton.");
- }
-}
More information about the CIG-COMMITS
mailing list