[cig-commits] commit 2448 by ian.rose to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Apr 9 12:31:26 PDT 2014
Revision 2448
Add the ability to do convection in a sphere, with associated gravity model
A trunk/aspect/cookbooks/future/sphere.prm
A trunk/aspect/include/aspect/geometry_model/sphere.h
U trunk/aspect/include/aspect/gravity_model/radial.h
A trunk/aspect/source/geometry_model/sphere.cc
U trunk/aspect/source/gravity_model/radial.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2448&peg=2448
Diff:
Added: trunk/aspect/cookbooks/future/sphere.prm
===================================================================
--- trunk/aspect/cookbooks/future/sphere.prm (rev 0)
+++ trunk/aspect/cookbooks/future/sphere.prm 2014-04-09 19:31:23 UTC (rev 2448)
@@ -0,0 +1,75 @@
+# A simple setup for convection in a sphere
+# manual for more information.
+set Dimension = 2
+set Use years in output instead of seconds = true
+set End time = 5.0e8
+set Output directory = output
+
+
+subsection Material model
+ set Model name = simple
+ subsection Simple model
+ set Thermal expansion coefficient = 4e-5
+ set Viscosity = 1e20
+ end
+end
+
+
+subsection Geometry model
+ set Model name = sphere
+ subsection Sphere
+ set Radius = 1250000
+ end
+end
+
+
+subsection Model settings
+ set Zero velocity boundary indicators =
+ set Tangential velocity boundary indicators = 0
+ set Prescribed velocity boundary indicators =
+ set Fixed temperature boundary indicators = 0
+ set Include shear heating = false
+ set Remove nullspace = net rotation
+end
+
+
+subsection Boundary temperature model
+ set Model name = constant
+ subsection Constant
+ set Boundary indicator to temperature mappings = 0 : 0.0
+ end
+end
+
+
+subsection Initial conditions
+ set Model name = function
+ subsection Function
+ set Function expression = if( sqrt( (x-5.e5)^2 + (y-5.e5)^2) < 1.e5, 0, 100)
+ end
+end
+
+
+subsection Gravity model
+ set Model name = radial linear
+ subsection Radial linear
+ set Magnitude at surface = 4.5
+ end
+end
+
+
+subsection Mesh refinement
+ set Initial global refinement = 5
+ set Initial adaptive refinement = 2
+ set Strategy = temperature
+ set Time steps between mesh refinement = 15
+end
+
+
+subsection Postprocess
+ set List of postprocessors = visualization, velocity statistics, temperature statistics, heat flux statistics
+ subsection Visualization
+ set Output format = vtu
+ set Time between graphical output = 1e6
+ set Number of grouped files = 0
+ end
+end
Added: trunk/aspect/include/aspect/geometry_model/sphere.h
===================================================================
--- trunk/aspect/include/aspect/geometry_model/sphere.h (rev 0)
+++ trunk/aspect/include/aspect/geometry_model/sphere.h 2014-04-09 19:31:23 UTC (rev 2448)
@@ -0,0 +1,97 @@
+/*
+ Copyright (C) 2014 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/>.
+*/
+#ifndef __aspect__geometry_model_sphere_h
+#define __aspect__geometry_model_sphere_h
+
+#include <aspect/geometry_model/interface.h>
+
+
+namespace aspect
+{
+ namespace GeometryModel
+ {
+ using namespace dealii;
+
+ template <int dim>
+ class Sphere : public Interface<dim>
+ {
+ public:
+ /**
+ * Generate a coarse mesh for the geometry described by this class.
+ */
+ virtual
+ void create_coarse_mesh (parallel::distributed::Triangulation<dim> &coarse_grid) const;
+
+ /**
+ * Return the typical length scale one would expect of features in this geometry,
+ * assuming realistic parameters.
+ *
+ * As discussed in the step-32 tutorial program, an appropriate length scale for
+ * this geometry is 10km, so we return $10^4$ here.
+ */
+ virtual
+ double length_scale () const;
+
+ virtual
+ double depth(const Point<dim> &position) const;
+
+ virtual
+ Point<dim> representative_point(const double depth) const;
+
+ virtual
+ double maximal_depth() const;
+
+ /**
+ * Return the set of boundary indicators that are used by this model. This
+ * information is used to determine what boundary indicators can be used in
+ * the input file.
+ *
+ * The sphere model only has one boundary, with indicator zero.
+ */
+ virtual
+ std::set<types::boundary_id>
+ get_used_boundary_indicators () const;
+
+ static
+ void
+ declare_parameters (ParameterHandler &prm);
+
+ virtual
+ void
+ parse_parameters (ParameterHandler &prm);
+
+ /**
+ * Return the radius of the sphere.
+ */
+ double
+ radius () const;
+
+ public:
+ /**
+ * Radius of the sphere
+ */
+ double R;
+
+ };
+ }
+}
+
+
+#endif
Modified: trunk/aspect/include/aspect/gravity_model/radial.h
===================================================================
--- trunk/aspect/include/aspect/gravity_model/radial.h 2014-04-09 14:59:11 UTC (rev 2447)
+++ trunk/aspect/include/aspect/gravity_model/radial.h 2014-04-09 19:31:23 UTC (rev 2448)
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2011, 2012 by the authors of the ASPECT code.
+ Copyright (C) 2014 by the authors of the ASPECT code.
This file is part of ASPECT.
@@ -23,6 +23,7 @@
#ifndef __aspect__gravity_model_radial_h
#define __aspect__gravity_model_radial_h
+#include <aspect/simulator.h>
#include <aspect/gravity_model/interface.h>
namespace aspect
@@ -91,6 +92,50 @@
virtual Tensor<1,dim> gravity_vector (const Point<dim> &position) const;
};
+ /**
+ * A class that describes gravity as a radial vector of linearly
+ * decreasing magnitude with depth. Meant for use in the Sphere
+ * geometry model, where you expect that kind of field.
+ *
+ * @ingroup GravityModels
+ */
+ template <int dim>
+ class RadialLinear : public Interface<dim>, public virtual SimulatorAccess<dim>
+ {
+ public:
+ /**
+ * Return the gravity vector as a function of position.
+ */
+ virtual Tensor<1,dim> gravity_vector (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);
+
+ /**
+ * Initialize the SimulatorAccess so this can have access
+ * to the geometry model
+ */
+ virtual void initialize(const Simulator<dim> &);
+
+ private:
+ /**
+ * Magnitude of the gravity vector at the surface, m/s^2
+ */
+ double magnitude_at_surface;
+
+ };
}
}
Added: trunk/aspect/source/geometry_model/sphere.cc
===================================================================
--- trunk/aspect/source/geometry_model/sphere.cc (rev 0)
+++ trunk/aspect/source/geometry_model/sphere.cc 2014-04-09 19:31:23 UTC (rev 2448)
@@ -0,0 +1,150 @@
+/*
+ Copyright (C) 2014 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/>.
+*/
+
+
+#include <aspect/geometry_model/sphere.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+
+
+namespace aspect
+{
+ namespace GeometryModel
+ {
+ template <int dim>
+ void
+ Sphere<dim>::
+ create_coarse_mesh (parallel::distributed::Triangulation<dim> &coarse_grid) const
+ {
+ GridGenerator::hyper_ball (coarse_grid,
+ Point<dim>(),
+ R);
+ static const HyperBallBoundary<dim> boundary_ball(Point<dim>(), R);
+ coarse_grid.set_boundary (0, boundary_ball);
+ }
+
+
+ template <int dim>
+ std::set<types::boundary_id>
+ Sphere<dim>::
+ get_used_boundary_indicators () const
+ {
+ const types::boundary_id s[] = { 0 };
+ return std::set<types::boundary_id>(&s[0],
+ &s[sizeof(s)/sizeof(s[0])]);
+ }
+
+
+ template <int dim>
+ double
+ Sphere<dim>::
+ length_scale () const
+ {
+ // as described in the first ASPECT paper, a length scale of
+ // 10km = 1e4m works well for the pressure scaling for earth
+ // sized spherical shells. use a length scale that
+ // yields this value for the R0,R1 corresponding to earth
+ // but otherwise scales like (R1-R0)
+ return 1e4 * maximal_depth() / (6336000.-3481000.);
+ }
+
+
+
+ template <int dim>
+ double
+ Sphere<dim>::depth(const Point<dim> &position) const
+ {
+ return std::min (std::max (R-position.norm(), 0.), maximal_depth());
+ }
+
+
+
+ template <int dim>
+ Point<dim>
+ Sphere<dim>::representative_point(const double depth) const
+ {
+ Point<dim> p;
+ p(dim-1) = std::min (std::max(R - depth, 0.), maximal_depth());
+ return p;
+ }
+
+
+
+ template <int dim>
+ double
+ Sphere<dim>::maximal_depth() const
+ {
+ return R;
+ }
+
+ template <int dim>
+ double Sphere<dim>::radius () const
+ {
+ return R;
+ }
+
+
+ template <int dim>
+ void
+ Sphere<dim>::declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Geometry model");
+ {
+ prm.enter_subsection("Sphere");
+ {
+ prm.declare_entry ("Radius", "6371000",
+ Patterns::Double (0),
+ "Radius of the sphere. Units: m.");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+
+
+ template <int dim>
+ void
+ Sphere<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Geometry model");
+ {
+ prm.enter_subsection("Sphere");
+ {
+ R = prm.get_double ("Radius");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace GeometryModel
+ {
+ ASPECT_REGISTER_GEOMETRY_MODEL(Sphere,
+ "sphere",
+ "Geometry model for sphere with a user specified radius.")
+ }
+}
Modified: trunk/aspect/source/gravity_model/radial.cc
===================================================================
--- trunk/aspect/source/gravity_model/radial.cc 2014-04-09 14:59:11 UTC (rev 2447)
+++ trunk/aspect/source/gravity_model/radial.cc 2014-04-09 19:31:23 UTC (rev 2448)
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2011, 2012 by the authors of the ASPECT code.
+ Copyright (C) 2014 by the authors of the ASPECT code.
This file is part of ASPECT.
@@ -17,7 +17,6 @@
along with ASPECT; see the file doc/COPYING. If not see
<http://www.gnu.org/licenses/>.
*/
-/* $Id$ */
#include <aspect/gravity_model/radial.h>
@@ -83,6 +82,60 @@
return -(1.245e-6 * r + 7.714e13/r/r) * p / r;
}
+
+// ----------------------------- RadialLinear ----------------------
+ template <int dim>
+ Tensor<1,dim>
+ RadialLinear<dim>::gravity_vector (const Point<dim> &p) const
+ {
+ if(p.norm() == 0.0) return Tensor<1,dim>();
+
+ double depth = this->get_geometry_model().depth(p);
+ Tensor<1,dim> grav = -magnitude_at_surface * p/p.norm() *
+ (1.0 - depth/this->get_geometry_model().maximal_depth());
+ return grav;
+ }
+
+ template <int dim>
+ void
+ RadialLinear<dim>::declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Gravity model");
+ {
+ prm.enter_subsection("Radial linear");
+ {
+ prm.declare_entry ("Magnitude at surface", "9.8",
+ Patterns::Double (0),
+ "Magnitude of the radial gravity vector"
+ "at the surface of the domain, m/s^2");
+ }
+ prm.leave_subsection ();
+ }
+ prm.leave_subsection ();
+ }
+
+
+ template <int dim>
+ void
+ RadialLinear<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Gravity model");
+ {
+ prm.enter_subsection("Radial linear");
+ {
+ magnitude_at_surface = prm.get_double ("Magnitude at surface");
+ }
+ prm.leave_subsection ();
+ }
+ prm.leave_subsection ();
+ }
+
+ //Need to do this in order to get the depth from the gravity model.
+ template <int dim>
+ void RadialLinear<dim>::initialize(const Simulator<dim> &simulation)
+ {
+ SimulatorAccess<dim>::initialize(simulation);
+ }
}
}
@@ -104,5 +157,11 @@
"the core-mantle boundary as well as at the surface and "
"in between is physically correct under the assumption "
"of a constant density.")
+
+ ASPECT_REGISTER_GRAVITY_MODEL(RadialLinear,
+ "radial linear",
+ "A gravity model which is radially inward, where the magnitude"
+ "decreases linearly with depth, as you would get with a constant"
+ "density spherical domain.")
}
}
More information about the CIG-COMMITS
mailing list