[cig-commits] r1283 - in branches/active_compositions: . cookbooks data doc/manual include/aspect/velocity_boundary_conditions source/velocity_boundary_conditions
dannberg at dealii.org
dannberg at dealii.org
Thu Oct 18 12:09:49 PDT 2012
Author: dannberg
Date: 2012-10-18 13:09:49 -0600 (Thu, 18 Oct 2012)
New Revision: 1283
Added:
branches/active_compositions/cookbooks/stokes.prm
branches/active_compositions/data/velocity-boundary-conditions/
branches/active_compositions/gplates.prm
branches/active_compositions/include/aspect/velocity_boundary_conditions/gplates.h
branches/active_compositions/source/velocity_boundary_conditions/gplates.cc
Modified:
branches/active_compositions/doc/manual/manual.tex
Log:
Merge from mainline and add cookbook stokes.prm, which also got a subsection in the manual.
Added: branches/active_compositions/cookbooks/stokes.prm
===================================================================
--- branches/active_compositions/cookbooks/stokes.prm (rev 0)
+++ branches/active_compositions/cookbooks/stokes.prm 2012-10-18 19:09:49 UTC (rev 1283)
@@ -0,0 +1,111 @@
+############### Global parameters
+
+set Dimension = 3
+
+set Start time = 0
+set End time = 2e13
+
+set Output directory = output
+set Use years in output instead of seconds = false
+
+############### Parameters describing the model
+
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 2890000
+ set Y extent = 2890000
+ set Z extent = 2890000
+ end
+end
+
+
+subsection Model settings
+ set Tangential velocity boundary indicators = 0,1,2,3,4,5
+ set Fixed temperature boundary indicators = 0,1
+ set Include shear heating = false
+end
+
+
+subsection Material model
+ set Model name = simple
+
+ subsection Simple model
+ set Reference density = 3300
+ set Reference temperature = 1
+ set Thermal expansion coefficient = 4e-5
+ set Viscosity = 1e22
+ end
+end
+
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 9.81
+ end
+end
+
+
+############### Parameters describing the temperature field
+
+subsection Boundary temperature model
+ set Model name = box
+
+ subsection Box
+ set Bottom temperature = 1
+ set Left temperature = 1
+ set Right temperature = 1
+ set Top temperature = 1
+ set Front temperature = 1
+ set Back temperature = 1
+ end
+end
+
+subsection Initial conditions
+ set Model name = function
+
+ subsection Function
+ set Function expression = 1
+ end
+end
+
+############### Parameters describing the compositional field
+
+subsection Compositional fields
+ set Number of fields = 1
+end
+
+subsection Compositional initial conditions
+ set Model name = function
+
+ subsection Function
+ set Variable names = x,y,z
+ set Function expression = if(sqrt((x-1445000)*(x-1445000)+(y-1445000)*(y-1445000)+(z-1445000)*(z-1445000)) > 200000,0,1)
+ end
+end
+
+############### Parameters describing the discretization
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 6
+ set Initial global refinement = 5
+ set Refinement fraction = 0.2
+ set Strategy = Velocity
+end
+
+
+############### Parameters describing the what to do with the solution
+
+subsection Postprocess
+ set List of postprocessors = visualization
+
+ subsection Visualization
+ set Number of grouped files = 0
+ set Output format = vtu
+ set Time between graphical output = 1e6
+ set List of output variables = density, viscosity
+ end
+end
Modified: branches/active_compositions/doc/manual/manual.tex
===================================================================
--- branches/active_compositions/doc/manual/manual.tex 2012-10-18 19:06:24 UTC (rev 1282)
+++ branches/active_compositions/doc/manual/manual.tex 2012-10-18 19:09:49 UTC (rev 1283)
@@ -3287,6 +3287,208 @@
end
\end{lstlisting}
+\subsubsection{The ``Stokes' law'' benchmark}
+\label{sec:benchmark-stokes_law}
+
+Stokes' law was derived by George Gabriel Stokes in 1851 and describes the frictional force
+a sphere experiences in a laminar flowing viscous medium. A setup for testing this
+law is a sphere with the radius r falling in a highly viscous fluid with lower density. Due to its
+higher density the sphere is accelerated by the gravitational force. While
+the frictional force increases with the velocity of the falling particle,
+the buoyancy force remains constant. Thus, at some time the forces will
+be balanced and the settling velocity of the sphere $v_s$ will remain constant:
+
+\begin{align}
+ \label{eq:stokes-law}
+ \underbrace{6 \pi \, \eta \, r \, v_s}_{\text{frictional force}} =
+ \underbrace{4/3 \pi \, r^3 \, \Delta\rho \, g}_{\text{buoyancy force}}
+\end{align}
+
+\begin{align*}
+ \eta \quad &\text{dynamic viscosity of the fluid}\\
+ \Delta\rho \quad &\text{density difference between sphere and fluid}\\
+ g \quad &\text{gravitational acceleration}
+\end{align*}
+
+The resulting settling velocity is given by:
+\begin{align}
+ \label{eq:stokes-velo}
+ v_s = \frac{2}{9} \frac{\Delta\rho \, r^2 \, g}{\eta}.
+\end{align}
+
+Because we do not take into account inertia in our numerical computation,
+the falling particle will reach the constant settling velocity right after
+the first timestep.
+For the setup of this benchmark, we chose the following parameters:
+\begin{align*}
+ \label{eq:stokes-parameters}
+ r &= 200 \, \text{km}\\
+ \Delta\rho &= 100 \, \text{kg}/\text{m}^3\\
+ \eta &= 10^{22} \, \text{Pa s}\\
+ g &= 9.81 \, \text{m}/\text{s}^2\\
+ \vspace{1mm}
+ \text{resulting } v_s &= 8.72 \cdot 10^{-10} \, \text{m}/\text{s}
+\end{align*}
+
+\begin{figure}
+ \begin{center}
+ \includegraphics[width=0.55\textwidth]{cookbooks/benchmarks/stokes-velocity}
+ \hfill
+ \includegraphics[width=0.44\textwidth]{cookbooks/benchmarks/stokes-density}
+ \caption{Stokes benchmark. Both figures show only a 2D slice of the
+ three-dimensional model.
+ Left: The compositional field and overlaid to it some velocity vectors.
+ The composition is 1 inside a sphere with the radius of 200 km and 0
+ outside of this sphere. As the velocity vectors show, the sphere sinks
+ in the viscous medium.
+ Right: The density distribution of the model. The compositional density
+ contrast of 100 kg$/\text{m}^3$ leads to a higher density inside of the
+ sphere.}
+ \label{fig:solcx}
+ \end{center}
+\end{figure}
+
+To run this benchmark, you can use the following input file:
+
+\begin{lstlisting}[frame=single,language=prmfile]
+############### Global parameters
+
+set Dimension = 3
+
+set Start time = 0
+set End time = 2e13
+
+set Output directory = output
+set Use years in output instead of seconds = false
+
+############### Parameters describing the model
+
+subsection Geometry model
+ set Model name = box
+
+ subsection Box
+ set X extent = 2890000
+ set Y extent = 2890000
+ set Z extent = 2890000
+ end
+end
+
+
+subsection Model settings
+ set Tangential velocity boundary indicators = 0,1,2,3,4,5
+ set Fixed temperature boundary indicators = 0,1
+ set Include shear heating = false
+end
+
+
+subsection Material model
+ set Model name = simple
+
+ subsection Simple model
+ set Reference density = 3300
+ set Reference temperature = 1
+ set Thermal expansion coefficient = 4e-5
+ set Viscosity = 1e22
+ end
+end
+
+
+subsection Gravity model
+ set Model name = vertical
+
+ subsection Vertical
+ set Magnitude = 9.81
+ end
+end
+
+
+############### Parameters describing the temperature field
+
+subsection Boundary temperature model
+ set Model name = box
+
+ subsection Box
+ set Bottom temperature = 1
+ set Left temperature = 1
+ set Right temperature = 1
+ set Top temperature = 1
+ set Front temperature = 1
+ set Back temperature = 1
+ end
+end
+
+subsection Initial conditions
+ set Model name = function
+
+ subsection Function
+ set Function expression = 1
+ end
+end
+
+############### Parameters describing the compositional field
+
+subsection Compositional fields
+ set Number of fields = 1
+end
+
+subsection Compositional initial conditions
+ set Model name = function
+
+ subsection Function
+ set Variable names = x,y,z
+ set Function expression = if(sqrt((x-1445000)*(x-1445000)+(y-1445000)*(y-1445000)
+ +(z-1445000)*(z-1445000)) > 200000,0,1)
+ end
+end
+
+############### Parameters describing the discretization
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 6
+ set Initial global refinement = 5
+ set Refinement fraction = 0.2
+ set Strategy = Velocity
+end
+
+
+############### Parameters describing the what to do with the solution
+
+subsection Postprocess
+ set List of postprocessors = visualization
+
+ subsection Visualization
+ set Number of grouped files = 0
+ set Output format = vtu
+ set Time between graphical output = 1e6
+ set List of output variables = density, viscosity
+ end
+end
+\end{lstlisting}
+
+\vspace{5mm}
+
+Finally, we want to know the settling velocity of the sphere in our numerical
+model. You can visualize the output in different ways, one of it being ParaView.
+Here you can calculate the average velocity of the sphere using the filters
+
+\begin{enumerate}
+ \item Threshold (Scalars: C\_1, Lower Threshold 0.5, Upper Threshold 1),
+ \item Integrate Variables,
+ \item Cell Data to Point Data,
+ \item Calculator (use the formula sqrt(velocity\_x\textasciicircum2+
+ velocity\_y\textasciicircum2+velocity\_z\textasciicircum2)/Volume).
+\end{enumerate}
+
+If you now look at
+the Calculator object in the Spreadsheet View, you can see the average sinking
+velocity of the sphere in the column ``Result'' and compare it to the theoretical
+value ($v_s = 8.72 \cdot 10^{-10} \, \text{m}/\text{s}$).
+In this case, this value should be about 8.865 $\cdot 10^{-10} \, \text{m}/\text{s}$.
+The difference between the analytical and the numerical solution can be explained by
+different points: In our case the sphere is viscous and not rigid, as in Stokes'
+law, and we have a finite box instead of an infinite medium.
+
+
\section{Extending \aspect}
\label{sec:extending}
Copied: branches/active_compositions/gplates.prm (from rev 1281, trunk/aspect/gplates.prm)
===================================================================
--- branches/active_compositions/gplates.prm (rev 0)
+++ branches/active_compositions/gplates.prm 2012-10-18 19:09:49 UTC (rev 1283)
@@ -0,0 +1,93 @@
+set Dimension = 3
+set End time = 4e14
+set Timing output frequency = 10
+set Use years in output instead of seconds = false
+
+subsection Geometry model
+ set Model name = spherical shell
+end
+
+subsection Gravity model
+ set Model name = radial constant
+
+ subsection Radial constant
+ set Magnitude = 9.81
+ end
+
+end
+
+subsection Material model
+ set Model name = simple
+
+ subsection Simple model
+ # The value of the constant viscosity. Units: $kg/m/s$.
+ set Viscosity = 1e22
+ end
+
+end
+
+
+subsection Mesh refinement
+ set Initial adaptive refinement = 0
+
+ set Initial global refinement = 1
+
+ set Refinement fraction = 0.0
+
+ set Coarsening fraction = 0.0
+end
+
+
+subsection Model settings
+ set Fixed temperature boundary indicators = 0,1
+
+ # A comma separated list of integers denoting those boundaries on which the
+ # velocity is tangential but prescribed, i.e., where external forces act to
+ # prescribe a particular velocity. This is often used to prescribe a
+ # velocity that equals that of overlying plates.
+ set Prescribed velocity boundary indicators = 1:gplates
+
+ #set Zero velocity boundary indicators = 1
+ set Tangential velocity boundary indicators = 0
+end
+
+subsection Boundary velocity model
+
+ subsection GPlates model
+ set Data directory = data/velocity-boundary-conditions/gplates/
+ set Velocity file name = time_dependent.%d.gpml
+ set Time step = 1e14
+ end
+
+end
+
+subsection Initial conditions
+ set Model name = function
+ subsection Function
+ set Function expression = 1600.0
+ end
+end
+
+subsection Boundary temperature model
+ set Model name = spherical constant
+
+ subsection Spherical constant
+ set Inner temperature = 4273
+ set Outer temperature = 273
+ end
+
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization,velocity statistics,temperature statistics,heat flux statistics, depth average
+
+ subsection Visualization
+ set Time between graphical output = 3.2e13
+ end
+
+ subsection Depth average
+ set Time between graphical output = 3.2e13
+ end
+end
+
+
Copied: branches/active_compositions/include/aspect/velocity_boundary_conditions/gplates.h (from rev 1281, trunk/aspect/include/aspect/velocity_boundary_conditions/gplates.h)
===================================================================
--- branches/active_compositions/include/aspect/velocity_boundary_conditions/gplates.h (rev 0)
+++ branches/active_compositions/include/aspect/velocity_boundary_conditions/gplates.h 2012-10-18 19:09:49 UTC (rev 1283)
@@ -0,0 +1,158 @@
+/*
+ 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: zero_velocity.h 1071 2012-08-01 16:50:02Z bangerth $ */
+
+
+#ifndef __aspect__velocity_boundary_conditions_gplates_h
+#define __aspect__velocity_boundary_conditions_gplates_h
+
+#include <aspect/velocity_boundary_conditions/interface.h>
+#include <aspect/simulator.h>
+
+
+namespace aspect
+{
+ namespace VelocityBoundaryConditions
+ {
+ using namespace dealii;
+
+ namespace internal
+ {
+ class GPlatesLookup;
+ }
+
+ /**
+ * A class that implements prescribed velocity boundary conditions
+ * determined from a GPlates input files.
+ *
+ * @ingroup VelocityBoundaryConditionsModels
+ */
+ template <int dim>
+ class GPlates : public Interface<dim>
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ GPlates ();
+
+ /**
+ * Return the boundary velocity as a function of position. For the
+ * current class, this function returns value from gplates.
+ */
+ Tensor<1,dim>
+ boundary_velocity (const Point<dim> &position) const;
+
+ /**
+ * Initialization function. This function is called once at the beginning
+ * of the program and loads the first set of GPlates files describing initial
+ * conditions at the start time.
+ */
+ void
+ initialize (const GeometryModel::Interface<dim> &geometry_model);
+
+ /**
+ * A function that is called at the beginning of each time step
+ * to indicate what the model time is for which the boundary
+ * values will next be evaluated. Also loads the next velocity
+ * files if necessary and outputs a warning if the end of the set of
+ * velocity files if reached.
+ */
+ void
+ set_current_time (const double time);
+
+ /**
+ * 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.
+ */
+ void
+ parse_parameters (ParameterHandler &prm);
+
+ private:
+ /**
+ * Pointer to the geometry object in use.
+ */
+ const GeometryModel::Interface<dim> *geometry_model;
+
+ /**
+ * A variable that stores the current time of the simulation. Derived
+ * classes can query this variable. It is set at the beginning of each
+ * time step.
+ */
+ double current_time;
+
+ /**
+ * A variable that stores the currently used velocity file of a series.
+ * It gets updated if necessary by set_current_time.
+ */
+ unsigned int current_time_step;
+
+ /**
+ * Directory in which the gplates velocity are present.
+ */
+ std::string data_directory;
+
+ /**
+ * First part of filename of velocity files. The files have to have the pattern
+ * velocity_file_name.n.gpml where n is the number of the current timestep (starts
+ * from 0).
+ */
+ std::string velocity_file_name;
+
+ /**
+ * Time in model units (depends on other model inputs) between two velocity files.
+ */
+ double time_step;
+
+ /**
+ * Weight between velocity file n and n+1 while the current time is between the
+ * two values t(n) and t(n+1).
+ */
+ double time_weight;
+
+ /**
+ * State whether we have time_dependent boundary conditions. Switched off after
+ * finding no more velocity files to suppress attempts to read in new files.
+ */
+ bool time_dependent;
+
+ /**
+ * Pointer to an object that reads and processes data we get from gplates files.
+ */
+ std_cxx1x::shared_ptr<internal::GPlatesLookup> lookup;
+
+ /**
+ * Create a filename out of the name template.
+ */
+ std::string
+ create_filename (const int timestep) const;
+ };
+ }
+}
+
+
+#endif
Copied: branches/active_compositions/source/velocity_boundary_conditions/gplates.cc (from rev 1281, trunk/aspect/source/velocity_boundary_conditions/gplates.cc)
===================================================================
--- branches/active_compositions/source/velocity_boundary_conditions/gplates.cc (rev 0)
+++ branches/active_compositions/source/velocity_boundary_conditions/gplates.cc 2012-10-18 19:09:49 UTC (rev 1283)
@@ -0,0 +1,476 @@
+/*
+ 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: zero_velocity.cc 1071 2012-08-01 16:50:02Z bangerth $ */
+
+
+#include <aspect/global.h>
+#include <aspect/velocity_boundary_conditions/gplates.h>
+#include <aspect/geometry_model/spherical_shell.h>
+
+#include <deal.II/base/parameter_handler.h>
+#include <deal.II/base/table.h>
+#include <fstream>
+#include <iostream>
+#include <boost/property_tree/xml_parser.hpp>
+#include <boost/property_tree/ptree.hpp>
+#include <boost/lexical_cast.hpp>
+
+
+namespace aspect
+{
+ namespace VelocityBoundaryConditions
+ {
+
+ namespace internal
+ {
+ /**
+ * GPlatesLookup handles all kinds of tasks around looking up a certain velocity boundary
+ * condition from a gplates .gpml file. This class keeps around the contents of two sets
+ * of files, corresponding to two instances in time where GPlates provides us with data;
+ * the boundary values at one particular time are interpolated between the two currently
+ * loaded data sets.
+ */
+ class GPlatesLookup
+ {
+ public:
+ /**
+ * Initialize all members and the two pointers referring to the actual velocities
+ */
+ GPlatesLookup()
+ {
+ velocity_vals.reinit(0,0);
+ old_velocity_vals.reinit(0,0);
+
+ velocity_values = &velocity_vals;
+ old_velocity_values = &old_velocity_vals;
+
+ delta_phi = 0.0;
+ delta_theta = 0.0;
+ }
+
+ /**
+ * Check whether a file named filename exists.
+ */
+ bool fexists(const std::string &filename) const
+ {
+ std::ifstream ifile(filename.c_str());
+ return ifile;
+ }
+
+ /**
+ * Loads a gplates .gpml velocity file. Throws an exception if
+ * the file does not exist.
+ */
+ void load_file(const std::string &filename)
+ {
+ using boost::property_tree::ptree;
+ ptree pt;
+
+ // Check whether file exists, we do not want to throw
+ // an exception in case it does not, because it could be by purpose
+ // (i.e. the end of the boundary condition is reached)
+ AssertThrow (fexists(filename),
+ ExcMessage (std::string("GPlates file <")
+ +
+ filename
+ +
+ "> not found!"));
+
+ // populate tree structure pt
+ read_xml(filename, pt);
+
+
+
+ const int n_points = pt.get_child("gpml:FeatureCollection.gml:featureMember.gpml:VelocityField.gml:domainSet.gml:MultiPoint").size();
+ const int n_phi = static_cast<int>(std::sqrt(2.*n_points));
+ const int n_theta = n_phi /2;
+
+ if ((delta_theta != 0.0) || (delta_phi != 0.0))
+ {
+ Assert((delta_theta - numbers::PI / (n_theta-1)) <= 1e-7, ExcMessage("Resolution changed during boundary conditions. Interpolation is not working correctly."));
+ Assert((delta_phi - 2*numbers::PI / n_phi) <= 1e-7, ExcMessage("Resolution changed during boundary conditions. Interpolation is not working correctly."));
+ }
+
+ delta_theta = numbers::PI / (n_theta-1);
+ delta_phi = 2*numbers::PI / n_phi;
+
+ // swap pointers to old and new values, we overwrite the old ones
+ // and the new ones become the old ones
+ std::swap( velocity_values, old_velocity_values);
+
+ (*velocity_values).reinit(n_theta,n_phi);
+
+ std::string velos = pt.get<std::string>("gpml:FeatureCollection.gml:featureMember.gpml:VelocityField.gml:rangeSet.gml:DataBlock.gml:tupleList");
+ std::stringstream in(velos, std::ios::in);
+ AssertThrow (in,
+ ExcMessage (std::string("Couldn't find velocities. Is file native gpml format for velocities?")));
+
+ unsigned int i = 0;
+ char sep;
+ while (!in.eof())
+ {
+ double vtheta, vphi;
+ const double cmyr_si = 3.1557e9;
+
+ in >> vtheta >> sep >> vphi;
+
+ if (in.eof())
+ break;
+
+ (*velocity_values)[i%n_theta][i/n_theta][0] = vtheta/cmyr_si;
+ (*velocity_values)[i%n_theta][i/n_theta][1] = vphi/cmyr_si;
+
+ i++;
+ }
+ }
+
+
+ /**
+ * Returns the computed surface velocity in cartesian coordinates. Takes
+ * as input the position and current time weight.
+ *
+ * @param position The current position to compute velocity @param time_weight A weighting between
+ * the two current timesteps n and n+1
+ */
+ template <int dim>
+ Tensor<1,dim> surface_velocity(const Point<dim> &position, const double time_weight) const
+ {
+ const Point<dim> scoord = spherical_surface_coordinates(position);
+
+ const double idtheta = get_idtheta(scoord(0));
+ const unsigned int idxtheta = static_cast<unsigned int>(idtheta);
+ const unsigned int nexttheta = idxtheta + 1;
+
+
+ const double idphi = get_idphi(scoord(1));
+ const unsigned int idxphi = static_cast<unsigned int>(idphi);
+ const unsigned int nextphi = (idxphi+1) % velocity_values->n_cols();
+
+
+ Assert(idxtheta<velocity_values->n_rows(), ExcMessage("not in range"));
+ Assert(idxphi <velocity_values->n_cols(), ExcMessage("not in range"));
+
+ Assert(nexttheta<velocity_values->n_rows(), ExcMessage("not in range"));
+ Assert(nextphi <velocity_values->n_cols(), ExcMessage("not in range"));
+
+ // compute the coordinates of this point in the
+ // reference cell between the data points
+ const double xi = idtheta-static_cast<double>(idxtheta);
+ const double eta = idphi-static_cast<double>(idxphi);
+
+ Assert ((0 <= xi) && (xi <= 1), ExcInternalError());
+ Assert ((0 <= eta) && (eta <= 1), ExcInternalError());
+ Assert ((0 <= time_weight) && (time_weight <= 1), ExcInternalError());
+
+ // use xi, eta and time_weight for a trilinear interpolation
+ // TODO: interpolation in cartesian probably more accurate
+
+
+ Tensor<1,2> surf_vel;
+
+ surf_vel[0] = time_weight *
+ ((1-xi)*(1-eta)*(*velocity_values)[idxtheta][idxphi][0] +
+ xi *(1-eta)*(*velocity_values)[nexttheta][idxphi][0] +
+ (1-xi)*eta *(*velocity_values)[idxtheta][nextphi][0] +
+ xi *eta *(*velocity_values)[nexttheta][nextphi][0]);
+
+ surf_vel[1] = time_weight *
+ ((1-xi)*(1-eta)*(*velocity_values)[idxtheta][idxphi][1] +
+ xi *(1-eta)*(*velocity_values)[nexttheta][idxphi][1] +
+ (1-xi)*eta *(*velocity_values)[idxtheta][nextphi][1] +
+ xi *eta *(*velocity_values)[nexttheta][nextphi][1]);
+
+
+ surf_vel[0] += (1-time_weight) *
+ ((1-xi)*(1-eta)*(*old_velocity_values)[idxtheta][idxphi][0] +
+ xi *(1-eta)*(*old_velocity_values)[nexttheta][idxphi][0] +
+ (1-xi)*eta *(*old_velocity_values)[idxtheta][nextphi][0] +
+ xi *eta *(*old_velocity_values)[nexttheta][nextphi][0]);
+
+ surf_vel[1] += (1-time_weight) *
+ ((1-xi)*(1-eta)*(*old_velocity_values)[idxtheta][idxphi][1] +
+ xi *(1-eta)*(*old_velocity_values)[nexttheta][idxphi][1] +
+ (1-xi)*eta *(*old_velocity_values)[idxtheta][nextphi][1] +
+ xi *eta *(*old_velocity_values)[nexttheta][nextphi][1]);
+
+
+ return sphere_to_cart_velocities(surf_vel,scoord);
+ }
+
+
+ private:
+
+ /**
+ * Tables which contain the velocities
+ */
+ dealii::Table<2,Tensor<1,2> > velocity_vals;
+ dealii::Table<2,Tensor<1,2> > old_velocity_vals;
+
+ /**
+ * Pointers to the actual tables.
+ * Used to avoid unnecessary copying
+ * of values.
+ */
+ dealii::Table<2,Tensor<1,2> > * velocity_values;
+ dealii::Table<2,Tensor<1,2> > * old_velocity_values;
+
+ /**
+ * Distances between adjacent point in the Lat/Lon grid
+ */
+ double delta_phi,delta_theta;
+
+ /**
+ * Returns spherical coordinates of a cartesian position.
+ */
+ template <int dim>
+ Point<dim> spherical_surface_coordinates(const Point<dim> &position) const
+ {
+ double coord[dim];
+ Point<dim> scoord;
+
+ if (dim == 3)
+ {
+ coord[0] = std::acos(position(2)/std::sqrt(position.square())); // Theta
+ coord[1] = std::atan2(position(1),position(0)); // Phi
+ if (coord[1] < 0.0) coord[1] = 2*numbers::PI + coord[1]; // correct phi to [0,2*pi]
+ coord[2] = std::sqrt(position.square()); // R
+
+ return Point<dim> (coord[0],coord[1],coord[2]);
+ }
+ else
+ return Point<dim>(); // TODO: define for 2 Dimensions
+ }
+
+ /**
+ * Returns cartesian velocities calculated from surface velocities and position in spherical coordinates
+ *
+ * @param s_velocities Surface velocities in spherical coordinates (theta, phi) @param s_position Position
+ * in spherical coordinates (theta,phi,radius)
+ */
+ template <int dim>
+ Tensor<1,dim> sphere_to_cart_velocities(const Tensor<1,2> &s_velocities, const Point<dim> &s_position) const
+ {
+ Tensor<1,dim> velocities;
+ if (dim == 3)
+ {
+ velocities[0] = -1.0 * s_velocities[1] * std::sin(s_position[1]) + s_velocities[0]*std::cos(s_position[0])*std::cos(s_position[1]);
+ velocities[1] = s_velocities[1]*std::cos(s_position[1])+s_velocities[0]*std::cos(s_position[0])*std::sin(s_position[1]);
+ velocities[2] = -1.0*s_velocities[0]*std::sin(s_position[0]);
+ return velocities;
+ }
+ else
+ return Tensor<1,dim>(); // TODO: define for 2 Dimensions
+ }
+
+ /**
+ * calculates the phi-index given a certain phi
+ */
+ double get_idphi(const double phi_) const
+ {
+ const double phi = std::max(std::min(phi_,2*numbers::PI-1e-7),0.0);
+
+ Assert(phi>=0, ExcMessage("not in range"));
+ Assert(phi<=2*numbers::PI, ExcMessage("not in range"));
+ return phi/delta_phi;
+ }
+
+ /**
+ * calculates the theta-index given a certain polar angle
+ */
+ double get_idtheta(const double theta_) const
+ {
+ double theta = std::max(std::min(theta_,numbers::PI-1e-7),0.0);
+
+ Assert(theta>=0, ExcMessage("not in range"));
+ Assert(theta<=numbers::PI, ExcMessage("not in range"));
+ return theta/delta_theta;
+ }
+
+
+ };
+ }
+
+
+ template <int dim>
+ GPlates<dim>::GPlates ()
+ :
+ lookup (new internal::GPlatesLookup)
+ {}
+
+
+
+ template <int dim>
+ void
+ GPlates<dim>::initialize (const GeometryModel::Interface<dim> &geometry_model)
+ {
+ int loaded;
+
+ Assert(dim==3, ExcMessage("gplates velocity model just works for dim = 3"));
+ Assert (dynamic_cast<const GeometryModel::SphericalShell<dim>*> (&geometry_model) != 0,
+ ExcMessage ("This boundary condition can only be used if the geometry "
+ "is a spherical shell."));
+
+ // Initialize variables
+ current_time_step = 0;
+ time_weight = 0.0;
+
+ // Load first velocity file
+ lookup->load_file(create_filename(current_time_step));
+
+ // Load next velocity file for interpolation
+ try
+ {
+ lookup->load_file(create_filename(current_time_step+1));
+ time_dependent = true;
+ }
+ catch (...)
+ // if loading the next time step's file failed, then simply
+ // take the last one a second time
+ {
+ // Next velocity file not found --> Constant velocities
+ lookup->load_file(create_filename(current_time_step));
+ time_dependent = false;
+
+ // Give warning if first processor
+ if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0)
+ std::cout << std::endl
+ << "Loading second velocity file did not succeed."
+ << "Assuming constant boundary conditions model." << std::endl
+ << std::endl;
+ }
+ }
+
+ template <int dim>
+ std::string
+ GPlates<dim>::create_filename (const int timestep) const
+ {
+ std::string templ = data_directory+velocity_file_name;
+
+ char filename[150];
+ snprintf (filename, 150, templ.c_str(), timestep);
+ return std::string (filename);
+ }
+
+ template <int dim>
+ void
+ GPlates<dim>::set_current_time (const double time)
+ {
+ current_time = time;
+ if (time_dependent)
+ {
+ Assert ((0 <= time_weight) && (time_weight <= 1),
+ ExcMessage("Error in set_current_time. Time_weight has to be in [0,1]"));
+
+ if (static_cast<unsigned int>(current_time / time_step) > current_time_step)
+ {
+ current_time_step = static_cast<unsigned int>(current_time / time_step);
+
+ // Load next velocity file for interpolation
+ try
+ {
+ lookup->load_file(create_filename(current_time_step+1));
+ }
+ catch (...)
+ // of loading failed, simply take the previous file a second time
+ {
+ // Next velocity file not found --> Constant velocities
+ lookup->load_file(create_filename(current_time_step));
+ // no longer consider the problem time dependent from here on out
+ time_dependent = false;
+
+ // Give warning if first processor
+ if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD)==0)
+ std::cout << std::endl << "Loading new velocity file did not succeed." << std::endl
+ << "Assuming constant boundary conditions for rest of model." << std::endl
+ << std::endl;
+ }
+ }
+ time_weight = current_time/time_step - current_time_step;
+
+ }
+ }
+
+
+ template <int dim>
+ Tensor<1,dim>
+ GPlates<dim>::
+ boundary_velocity (const Point<dim> &position) const
+ {
+ return lookup->surface_velocity(position,time_weight);
+ }
+
+
+ template <int dim>
+ void
+ GPlates<dim>::declare_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Boundary velocity model");
+ {
+ prm.enter_subsection("GPlates model");
+ {
+ prm.declare_entry ("Data directory", "data/velocity-boundary-conditions/gplates/",
+ Patterns::DirectoryName (),
+ "The path to the model data.");
+ prm.declare_entry ("Velocity file name", "phi_test.%d",
+ Patterns::Anything (),
+ "The file name of the material data. Provide file in format: "
+ "(Velocity file name).%d.gpml where %d is any sprintf integer qualifier, "
+ "specifying the format of the current file number.");
+ prm.declare_entry ("Time step", "3.1558e13",
+ Patterns::Anything (),
+ "Time step between following velocity files."
+ " Default is one million years expressed in SI units.");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+ template <int dim>
+ void
+ GPlates<dim>::parse_parameters (ParameterHandler &prm)
+ {
+ prm.enter_subsection("Boundary velocity model");
+ {
+ prm.enter_subsection("GPlates model");
+ {
+ data_directory = prm.get ("Data directory");
+ velocity_file_name = prm.get ("Velocity file name");
+ time_step = prm.get_double ("Time step");
+ }
+ prm.leave_subsection();
+ }
+ prm.leave_subsection();
+ }
+
+ }
+}
+
+// explicit instantiations
+namespace aspect
+{
+ namespace VelocityBoundaryConditions
+ {
+ ASPECT_REGISTER_VELOCITY_BOUNDARY_CONDITIONS(GPlates,
+ "gplates",
+ "Implementation of a model in which the boundary "
+ "velocity is derived from GPlates.");
+ }
+}
More information about the CIG-COMMITS
mailing list