[cig-commits] commit 2314 by gassmoeller to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Thu Feb 20 12:28:02 PST 2014
Revision 2314
Included an initial condition plugin for temperature that produces a harmonic perturbation of a constant/adiabatic (incompressible/compressible) temperature profile for arbitrary geometries and dimensions. It uses spherical harmonics for 3D spherical shells.
A trunk/aspect/include/aspect/initial_conditions/harmonic_perturbation.h
A trunk/aspect/source/initial_conditions/harmonic_perturbation.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2314&peg=2314
Diff:
Added: trunk/aspect/include/aspect/initial_conditions/harmonic_perturbation.h
===================================================================
--- trunk/aspect/include/aspect/initial_conditions/harmonic_perturbation.h (rev 0)
+++ trunk/aspect/include/aspect/initial_conditions/harmonic_perturbation.h 2014-02-20 20:27:59 UTC (rev 2314)
@@ -0,0 +1,116 @@
+/*
+ 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: spherical_shell.h 853 2012-03-24 21:52:36Z bangerth $ */
+
+
+#ifndef __aspect__initial_conditions_harmonic_perturbation_h
+#define __aspect__initial_conditions_harmonic_perturbation_h
+
+#include <aspect/initial_conditions/interface.h>
+#include <aspect/simulator.h>
+
+
+namespace aspect
+{
+ namespace InitialConditions
+ {
+ using namespace dealii;
+
+ /**
+ * A class that describes a perturbed initially constant temperature field for
+ * any geometry model or dimension in shape of a harmonic function.
+ * For 3D spherical shell models this is achieved by using spherical harmonics,
+ * in any other case sine function are scaled to fit the model geometry.
+ *
+ * @ingroup InitialConditionsModels
+ */
+ template <int dim>
+ class HarmonicPerturbation : public Interface<dim>, public ::aspect::SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * Return the initial temperature as a function of position.
+ */
+ virtual
+ double initial_temperature (const Point<dim> &position) const;
+
+ /**
+ * Declare the parameters this class takes through input files.
+ */
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ /**
+ * Read the parameters this class declares from the parameter
+ * file.
+ */
+ virtual
+ void
+ parse_parameters (ParameterHandler &prm);
+
+
+ private:
+
+ /**
+ * Returns spherical coordinates of a cartesian position.
+ */
+ const Tensor<1,dim>
+ spherical_surface_coordinates(const Tensor<1,dim> &position) const;
+
+ /**
+ * The radial/depth wave number of the harmonic perturbation.
+ * All wave number variables are in fact twice the wave number in a
+ * mathematical sense. This allows the user to prescribe
+ * a single up- / downswing or half periods.
+ */
+ int vertical_wave_number;
+
+ /**
+ * 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 maximal magnitude of the harmonic perturbation.
+ */
+ double magnitude;
+
+ /**
+ * The background temperature the harmonic perturbation is applied on in
+ * an incompressible material model. In case of a compressible material model
+ * the perturbation is applied on top of an adiabatic profile and this variable
+ * is not used at all.
+ */
+ double reference_temperature;
+ };
+ }
+}
+
+#endif
Added: trunk/aspect/source/initial_conditions/harmonic_perturbation.cc
===================================================================
--- trunk/aspect/source/initial_conditions/harmonic_perturbation.cc (rev 0)
+++ trunk/aspect/source/initial_conditions/harmonic_perturbation.cc 2014-02-20 20:27:59 UTC (rev 2314)
@@ -0,0 +1,203 @@
+/*
+ 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: spherical_shell.cc 1651 2013-04-29 00:20:56Z heister $ */
+
+
+#include <aspect/initial_conditions/harmonic_perturbation.h>
+#include <aspect/geometry_model/box.h>
+#include <aspect/geometry_model/spherical_shell.h>
+
+#include <boost/math/special_functions/spherical_harmonic.hpp>
+
+namespace aspect
+{
+ namespace InitialConditions
+ {
+ template <int dim>
+ double
+ HarmonicPerturbation<dim>::
+ initial_temperature (const Point<dim> &position) const
+ {
+
+ // use either the user-input reference temperature as background temperature
+ // (incompressible model) or the adiabatic temperature profile (compressible model)
+ const double background_temperature = this->get_material_model().is_compressible() ?
+ this->adiabatic_conditions->temperature(position) :
+ reference_temperature;
+
+ // s = fraction of the way from
+ // the inner to the outer
+ // boundary; 0<=s<=1
+ const double s = this->geometry_model->depth(position) / this->geometry_model->maximal_depth();
+
+ const double depth_perturbation = std::sin(vertical_wave_number*s*numbers::PI);
+
+
+ double lateral_perturbation = 0.0;
+
+ if (const GeometryModel::SphericalShell<dim> *
+ geometry_model = dynamic_cast <const GeometryModel::SphericalShell<dim>*> (this->geometry_model))
+ {
+ // 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 = 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]);
+ }
+ }
+ else if (const GeometryModel::Box<dim> *
+ geometry_model = dynamic_cast <const GeometryModel::Box<dim>*> (this->geometry_model))
+ {
+ // In case of Box model use a sine as lateral perturbation
+ // that is scaled to the extent of the geometry.
+ // This way the perturbation is alway 0 at the model borders.
+ const Point<dim> extent = geometry_model->get_extents();
+
+ if (dim==2)
+ {
+ lateral_perturbation = std::sin(lateral_wave_number_1*position(0)*numbers::PI/extent(0));
+ }
+ else if (dim==3)
+ {
+ lateral_perturbation = std::sin(lateral_wave_number_1*position(0)*numbers::PI/extent(0))
+ * std::sin(lateral_wave_number_2*position(1)*numbers::PI/extent(1));
+ }
+ }
+ else
+ AssertThrow (false,
+ ExcMessage ("Not a valid geometry model for the initial conditions model"
+ "harmonic perturbation."));
+
+ return background_temperature + magnitude * depth_perturbation * lateral_perturbation;
+ }
+
+ template <int dim>
+ const Tensor<1,dim>
+ HarmonicPerturbation<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
+ HarmonicPerturbation<dim>::declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Initial conditions");
+ {
+ prm.enter_subsection("Harmonic perturbation");
+ {
+ prm.declare_entry ("Vertical wave number", "1",
+ Patterns::Integer (),
+ "Doubled radial wave number of the harmonic perturbation. "
+ " 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.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.declare_entry ("Magnitude", "1.0",
+ Patterns::Double (0),
+ "The magnitude of the Harmonic perturbation.");
+ prm.declare_entry ("Reference temperature", "1600.0",
+ Patterns::Double (0),
+ "The reference temperature that is perturbed by the"
+ "harmonic function. Only used in incompressible models.");
+ }
+ prm.leave_subsection ();
+ }
+ prm.leave_subsection ();
+ }
+
+
+ template <int dim>
+ void
+ HarmonicPerturbation<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Initial conditions");
+ {
+ prm.enter_subsection("Harmonic perturbation");
+ {
+ vertical_wave_number = prm.get_integer ("Vertical wave number");
+ lateral_wave_number_1 = prm.get_integer ("Lateral wave number one");
+ lateral_wave_number_2 = prm.get_integer ("Lateral wave number two");
+ magnitude = prm.get_double ("Magnitude");
+ reference_temperature = prm.get_double ("Reference temperature");
+ }
+ prm.leave_subsection ();
+ }
+ prm.leave_subsection ();
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace InitialConditions
+ {
+ ASPECT_REGISTER_INITIAL_CONDITIONS(HarmonicPerturbation,
+ "harmonic perturbation",
+ "An initial temperature field in which the temperature "
+ "is perturbed following a harmonic function (spherical "
+ "harmonic or sine depending on geometry and dimension) "
+ "in lateral and radial direction from an otherwise "
+ "constant temperature (incompressible model) or adiabatic "
+ "reference profile (compressible model).")
+ }
+}
More information about the CIG-COMMITS
mailing list