[cig-commits] r1297 - in trunk/aspect: . cookbooks doc/manual doc/manual/cookbooks/benchmarks include/aspect include/aspect/compositional_initial_conditions include/aspect/material_model source source/compositional_initial_conditions source/material_model source/postprocess source/postprocess/visualization source/simulator
dannberg at dealii.org
dannberg at dealii.org
Fri Oct 19 08:30:59 PDT 2012
Author: dannberg
Date: 2012-10-19 09:30:58 -0600 (Fri, 19 Oct 2012)
New Revision: 1297
Added:
trunk/aspect/cookbooks/stokes.prm
trunk/aspect/doc/manual/cookbooks/benchmarks/stokes-density.png
trunk/aspect/doc/manual/cookbooks/benchmarks/stokes-velocity.png
Modified:
trunk/aspect/
trunk/aspect/doc/manual/manual.tex
trunk/aspect/include/aspect/adiabatic_conditions.h
trunk/aspect/include/aspect/compositional_initial_conditions/interface.h
trunk/aspect/include/aspect/material_model/duretz_et_al.h
trunk/aspect/include/aspect/material_model/interface.h
trunk/aspect/include/aspect/material_model/simple.h
trunk/aspect/include/aspect/material_model/steinberger.h
trunk/aspect/include/aspect/material_model/table.h
trunk/aspect/include/aspect/material_model/tan_gurnis.h
trunk/aspect/include/aspect/simulator.h
trunk/aspect/source/adiabatic_conditions.cc
trunk/aspect/source/compositional_initial_conditions/
trunk/aspect/source/compositional_initial_conditions/interface.cc
trunk/aspect/source/material_model/duretz_et_al.cc
trunk/aspect/source/material_model/interface.cc
trunk/aspect/source/material_model/simple.cc
trunk/aspect/source/material_model/steinberger.cc
trunk/aspect/source/material_model/table.cc
trunk/aspect/source/material_model/tan_gurnis.cc
trunk/aspect/source/postprocess/depth_average.cc
trunk/aspect/source/postprocess/heat_flux_statistics.cc
trunk/aspect/source/postprocess/table_heat_flux_statistics.cc
trunk/aspect/source/postprocess/table_velocity_statistics.cc
trunk/aspect/source/postprocess/velocity_statistics.cc
trunk/aspect/source/postprocess/visualization.cc
trunk/aspect/source/postprocess/visualization/density.cc
trunk/aspect/source/postprocess/visualization/friction_heating.cc
trunk/aspect/source/postprocess/visualization/nonadiabatic_pressure.cc
trunk/aspect/source/postprocess/visualization/nonadiabatic_temperature.cc
trunk/aspect/source/postprocess/visualization/seismic_vp.cc
trunk/aspect/source/postprocess/visualization/seismic_vs.cc
trunk/aspect/source/postprocess/visualization/specific_heat.cc
trunk/aspect/source/postprocess/visualization/strain_rate.cc
trunk/aspect/source/postprocess/visualization/thermal_expansivity.cc
trunk/aspect/source/postprocess/visualization/thermodynamic_phase.cc
trunk/aspect/source/postprocess/visualization/viscosity.cc
trunk/aspect/source/postprocess/visualization/viscosity_ratio.cc
trunk/aspect/source/simulator/assembly.cc
trunk/aspect/source/simulator/core.cc
trunk/aspect/source/simulator/helper_functions.cc
trunk/aspect/source/simulator/parameters.cc
trunk/aspect/source/simulator/simulator_access.cc
Log:
Reintegrate branch active_compositions: material parameters in dependence of the compositional fields
Property changes on: trunk/aspect
___________________________________________________________________
Modified: svn:mergeinfo
- /branches/compositional:1141-1251
/branches/fully-nonlinear:542-728
+ /branches/active_compositions:1257-1296
/branches/compositional:1141-1251
/branches/fully-nonlinear:542-728
Copied: trunk/aspect/cookbooks/stokes.prm (from rev 1296, branches/active_compositions/cookbooks/stokes.prm)
===================================================================
--- trunk/aspect/cookbooks/stokes.prm (rev 0)
+++ trunk/aspect/cookbooks/stokes.prm 2012-10-19 15:30:58 UTC (rev 1297)
@@ -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
Copied: trunk/aspect/doc/manual/cookbooks/benchmarks/stokes-density.png (from rev 1296, branches/active_compositions/doc/manual/cookbooks/benchmarks/stokes-density.png)
===================================================================
(Binary files differ)
Copied: trunk/aspect/doc/manual/cookbooks/benchmarks/stokes-velocity.png (from rev 1296, branches/active_compositions/doc/manual/cookbooks/benchmarks/stokes-velocity.png)
===================================================================
(Binary files differ)
Modified: trunk/aspect/doc/manual/manual.tex
===================================================================
--- trunk/aspect/doc/manual/manual.tex 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/doc/manual/manual.tex 2012-10-19 15:30:58 UTC (rev 1297)
@@ -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}
Modified: trunk/aspect/include/aspect/adiabatic_conditions.h
===================================================================
--- trunk/aspect/include/aspect/adiabatic_conditions.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/adiabatic_conditions.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -27,6 +27,7 @@
#include <aspect/material_model/interface.h>
#include <aspect/geometry_model/interface.h>
#include <aspect/gravity_model/interface.h>
+#include <aspect/compositional_initial_conditions/interface.h>
#include <deal.II/base/point.h>
@@ -56,8 +57,10 @@
AdiabaticConditions (const GeometryModel::Interface<dim> &geometry_model,
const GravityModel::Interface<dim> &gravity_model,
const MaterialModel::Interface<dim> &material_model,
+ const CompositionalInitialConditions::Interface<dim> &compositional_initial_conditions,
const double surface_pressure,
- const double surface_temperature);
+ const double surface_temperature,
+ const unsigned int n_compositional_fields);
/**
* Return the adiabatic temperature at a given point of the domain.
Modified: trunk/aspect/include/aspect/compositional_initial_conditions/interface.h
===================================================================
--- trunk/aspect/include/aspect/compositional_initial_conditions/interface.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/compositional_initial_conditions/interface.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -26,7 +26,6 @@
#include <aspect/plugins.h>
#include <aspect/geometry_model/interface.h>
#include <aspect/boundary_temperature/interface.h>
-#include <aspect/adiabatic_conditions.h>
#include <deal.II/base/point.h>
#include <deal.II/base/parameter_handler.h>
@@ -64,9 +63,7 @@
* adiabatic conditions and store them so that derived classes can access them.
*/
void
- initialize (const GeometryModel::Interface<dim> &geometry_model,
- const BoundaryTemperature::Interface<dim> &boundary_temperature,
- const AdiabaticConditions<dim> &adiabatic_conditions);
+ initialize (const GeometryModel::Interface<dim> &geometry_model);
/**
* Return the initial temperature as a function of position.
@@ -101,16 +98,6 @@
*/
const GeometryModel::Interface<dim> *geometry_model;
- /**
- * Pointer to an object that described the boundary values
- * for the temperature field.
- */
- const BoundaryTemperature::Interface<dim> *boundary_temperature;
-
- /**
- * Pointer to an object that describes adiabatic conditions.
- */
- const AdiabaticConditions<dim> *adiabatic_conditions;
};
@@ -150,9 +137,7 @@
template <int dim>
Interface<dim> *
create_initial_conditions (ParameterHandler &prm,
- const GeometryModel::Interface<dim> &geometry_model,
- const BoundaryTemperature::Interface<dim> &boundary_temperature,
- const AdiabaticConditions<dim> &adiabatic_conditions);
+ const GeometryModel::Interface<dim> &geometry_model);
/**
Modified: trunk/aspect/include/aspect/material_model/duretz_et_al.h
===================================================================
--- trunk/aspect/include/aspect/material_model/duretz_et_al.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/material_model/duretz_et_al.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -76,27 +76,33 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_conductivity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
* @}
@@ -244,27 +250,33 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_conductivity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
* @}
@@ -365,27 +377,33 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_conductivity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
* @}
Modified: trunk/aspect/include/aspect/material_model/interface.h
===================================================================
--- trunk/aspect/include/aspect/material_model/interface.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/material_model/interface.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -86,7 +86,7 @@
*/
/**
* Return the viscosity $\eta$ of the model as a function of temperature,
- * pressure, strain rate, and position.
+ * pressure, composition, strain rate, and position.
*
* @note The strain rate given as the third argument of this function
* is computed as $\varepsilon(\mathbf u)=\frac 12 (\nabla \mathbf u +
@@ -98,6 +98,7 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const = 0;
@@ -107,6 +108,7 @@
*/
virtual double viscosity_ratio (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strainrate,
const Point<dim> &position) const;
@@ -116,6 +118,7 @@
*/
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const = 0;
/**
@@ -125,6 +128,7 @@
*/
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const = 0;
/**
@@ -133,6 +137,7 @@
*/
virtual double specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const = 0;
/**
@@ -149,6 +154,7 @@
*/
virtual double thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
@@ -167,6 +173,7 @@
*/
virtual double thermal_conductivity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const = 0;
/**
* @}
@@ -219,7 +226,7 @@
* Return whether the model is compressible or not. Incompressibility
* does not necessarily imply that the density is constant; rather, it
* may still depend on temperature or pressure. In the current
- * context, compressibility means whether we should solve the contuity
+ * context, compressibility means whether we should solve the continuity
* equation as $\nabla \cdot (\rho \mathbf u)=0$ (compressible Stokes)
* or as $\nabla \cdot \mathbf{u}=0$ (incompressible Stokes).
*/
@@ -245,6 +252,7 @@
virtual double
viscosity_derivative (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position,
const NonlinearDependence::Dependence dependence) const;
@@ -259,6 +267,7 @@
virtual double
density_derivative (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position,
const NonlinearDependence::Dependence dependence) const;
@@ -273,6 +282,7 @@
virtual double
compressibility_derivative (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position,
const NonlinearDependence::Dependence dependence) const;
@@ -287,6 +297,7 @@
virtual double
specific_heat_derivative (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position,
const NonlinearDependence::Dependence dependence) const;
@@ -301,6 +312,7 @@
virtual double
thermal_conductivity_derivative (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position,
const NonlinearDependence::Dependence dependence) const;
/**
@@ -360,6 +372,7 @@
double
seismic_Vp (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
@@ -377,6 +390,7 @@
double
seismic_Vs (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
* Return the Phase number of the model as a function of
@@ -391,11 +405,38 @@
virtual
unsigned int
thermodynamic_phase (const double temperature,
- const double pressure) const;
+ const double pressure,
+ const std::vector<double> &compositional_fields) const;
/**
* @}
*/
+ struct MaterialModelInputs
+ {
+ MaterialModelInputs(unsigned int n_points, unsigned int n_comp);
+
+ std::vector<Point<dim> > position;
+ std::vector<double> temperature;
+ std::vector<double> pressure;
+ std::vector<std::vector<double> > composition;
+ std::vector<SymmetricTensor<2,dim> > strain_rate;
+ };
+
+ struct MaterialModelOutputs
+ {
+ MaterialModelOutputs(unsigned int n_points);
+
+ std::vector<double> viscosities;
+ std::vector<double> densities;
+ std::vector<double> thermal_expansion_coefficients;
+ std::vector<double> specific_heat;
+ std::vector<double> thermal_conductivities;
+ std::vector<double> compressibilities;
+ bool is_compressible;
+ };
+
+ virtual void compute_parameters(const MaterialModelInputs & in, MaterialModelOutputs & out);
+
/**
* @name Functions used in dealing with run-time parameters
* @{
Modified: trunk/aspect/include/aspect/material_model/simple.h
===================================================================
--- trunk/aspect/include/aspect/material_model/simple.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/material_model/simple.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -24,6 +24,7 @@
#define __aspect__model_simple_h
#include <aspect/material_model/interface.h>
+#include <aspect/simulator.h>
namespace aspect
{
@@ -43,7 +44,7 @@
* @ingroup MaterialModels
*/
template <int dim>
- class Simple : public MaterialModel::Interface<dim>
+ class Simple : public MaterialModel::Interface<dim>, public ::aspect::SimulatorAccess<dim>
{
public:
/**
@@ -52,27 +53,33 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_conductivity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
* @}
Modified: trunk/aspect/include/aspect/material_model/steinberger.h
===================================================================
--- trunk/aspect/include/aspect/material_model/steinberger.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/material_model/steinberger.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -56,35 +56,43 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_conductivity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double seismic_Vp (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double seismic_Vs (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
* @}
Modified: trunk/aspect/include/aspect/material_model/table.h
===================================================================
--- trunk/aspect/include/aspect/material_model/table.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/material_model/table.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -48,28 +48,34 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double viscosity_ratio (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_conductivity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
* @}
@@ -180,6 +186,7 @@
*/
virtual double thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
@@ -207,7 +214,8 @@
* quantities.
*/
virtual double seismic_Vp (const double temperature,
- const double pressure) const;
+ const double pressure,
+ const std::vector<double> &compositional_fields) const;
/**
* the seismic shear wave speed
@@ -217,7 +225,8 @@
* quantities.
*/
virtual double seismic_Vs (const double temperature,
- const double pressure) const;
+ const double pressure,
+ const std::vector<double> &compositional_fields) const;
/**
* the phase of the composition at given pressure and temperature
@@ -229,7 +238,8 @@
* quantities.
*/
virtual unsigned int thermodynamic_phase (const double temperature,
- const double pressure) const;
+ const double pressure,
+ const std::vector<double> &compositional_fields) const;
/**
Modified: trunk/aspect/include/aspect/material_model/tan_gurnis.h
===================================================================
--- trunk/aspect/include/aspect/material_model/tan_gurnis.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/material_model/tan_gurnis.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -51,27 +51,33 @@
*/
virtual double viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const;
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
virtual double thermal_conductivity (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
* @}
Modified: trunk/aspect/include/aspect/simulator.h
===================================================================
--- trunk/aspect/include/aspect/simulator.h 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/include/aspect/simulator.h 2012-10-19 15:30:58 UTC (rev 1297)
@@ -362,6 +362,19 @@
const BoundaryTemperature::Interface<dim> &
get_boundary_temperature () const;
+ /**
+ * Copy the values of the compositional fields at the quadrature point
+ * q given as input parameter to the output vector
+ * composition_values_at_q_point.
+ *
+ * This function is implemented in
+ * <code>source/simulator/helper_functions.cc</code>.
+ */
+ void
+ get_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
+ const unsigned int q,
+ std::vector<double> &composition_values_at_q_point) const;
+
/** @} */
private:
@@ -1161,6 +1174,11 @@
* given the values and gradients of the solution passed as
* arguments.
*
+ * @param index The index of the field we need the artificial
+ * diffusion coefficient for:
+ * 0 temperature
+ * 1...n_compositional_fields compositional field
+ *
* This function is implemented in
* <code>source/simulator/assembly.cc</code>.
*/
@@ -1177,14 +1195,107 @@
const std::vector<SymmetricTensor<2,dim> > &old_old_strain_rates,
const std::vector<double> &old_pressure,
const std::vector<double> &old_old_pressure,
+ const std::vector<std::vector<double>> &old_composition,
+ const std::vector<std::vector<double>> &old_old_composition,
const double global_u_infty,
const double global_T_variation,
const double average_temperature,
const double global_entropy_variation,
const std::vector<Point<dim> > &evaluation_points,
- const double cell_diameter) const;
+ const double cell_diameter,
+ const unsigned int index) const;
/**
+ * Compute the residual of the temperature equation to be used
+ * for the artificial diffusion coefficient value on a cell
+ * given the values and gradients of the solution
+ * passed as arguments.
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ void
+ compute_temperature_system_residual(const std::vector<double> &old_temperature,
+ const std::vector<double> &old_old_temperature,
+ const std::vector<Tensor<1,dim> > &old_temperature_grads,
+ const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
+ const std::vector<double> &old_temperature_laplacians,
+ const std::vector<double> &old_old_temperature_laplacians,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
+ const std::vector<SymmetricTensor<2,dim> > &old_strain_rates,
+ const std::vector<SymmetricTensor<2,dim> > &old_old_strain_rates,
+ const std::vector<double> &old_pressure,
+ const std::vector<double> &old_old_pressure,
+ const std::vector<std::vector<double>> &old_composition,
+ const std::vector<std::vector<double>> &old_old_composition,
+ const double average_temperature,
+ const std::vector<Point<dim> > &evaluation_points,
+ double &max_residual,
+ double &max_velocity) const;
+
+ /**
+ * Compute the residual of the composition equation to be used
+ * for the artificial diffusion coefficient value on a cell
+ * given the values and gradients of the solution
+ * passed as arguments.
+ *
+ * @param composition_index The index of the compositional field whose
+ * residual we want to get (0 <= composition_index < number of
+ * compositional fields in this problem).
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ void
+ compute_composition_system_residual(const std::vector<double> &old_composition,
+ const std::vector<double> &old_old_composition,
+ const std::vector<Tensor<1,dim> > &old_composition_grads,
+ const std::vector<Tensor<1,dim> > &old_old_composition_grads,
+ const std::vector<double> &old_composition_laplacians,
+ const std::vector<double> &old_old_composition_laplacians,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
+ const double average_composition,
+ const unsigned int composition_index,
+ double &max_residual,
+ double &max_velocity) const;
+ /**
+ * Extract the values of temperature, pressure, composition and
+ * optional strain rate for the current linearization point.
+ * These values are stored as input arguments for the material
+ * model. The compositional fields are extracted with the
+ * individual compositional fields as outer vectors and the values
+ * at each quadrature point as inner vectors, but the material
+ * model needs it the other way round. Hence, this vector of vectors
+ * is transposed.
+ *
+ * @param compute_strainrate If the strain rate should be computed
+ * or not
+ *
+ * This function is implemented in
+ * <code>source/simulator/assembly.cc</code>.
+ */
+ void
+ compute_material_model_input_values (const TrilinosWrappers::MPI::BlockVector ¤t_linearization_point,
+ const FEValues<dim,dim> &input_finite_element_values,
+ const bool compute_strainrate,
+ typename MaterialModel::Interface<dim>::MaterialModelInputs &material_model_inputs) const;
+
+ /**
+ * Copy the values of the compositional fields at the quadrature point
+ * q given as input parameter to the output vector
+ * composition_values_at_q_point.
+ *
+ * This function is implemented in
+ * <code>source/simulator/helper_functions.cc</code>.
+ */
+ void
+ extract_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
+ const unsigned int q,
+ std::vector<double> &composition_values_at_q_point) const;
+
+ /**
* Generate and output some statistics like timing information and memory consumption.
* Whether this function does anything or not is controlled through the
* variable aspect::output_parallel_statistics.
@@ -1253,7 +1364,7 @@
const std::auto_ptr<GravityModel::Interface<dim> > gravity_model;
const std::auto_ptr<BoundaryTemperature::Interface<dim> > boundary_temperature;
std::auto_ptr< InitialConditions::Interface<dim> > initial_conditions;
- std::auto_ptr<const CompositionalInitialConditions::Interface<dim> > compositional_initial_conditions;
+ std::auto_ptr<CompositionalInitialConditions::Interface<dim> > compositional_initial_conditions;
std::auto_ptr<const AdiabaticConditions<dim> > adiabatic_conditions;
std::map<types::boundary_id_t,std_cxx1x::shared_ptr<VelocityBoundaryConditions::Interface<dim> > > velocity_boundary_conditions;
/**
Modified: trunk/aspect/source/adiabatic_conditions.cc
===================================================================
--- trunk/aspect/source/adiabatic_conditions.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/adiabatic_conditions.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -34,8 +34,10 @@
AdiabaticConditions<dim>::AdiabaticConditions(const GeometryModel::Interface<dim> &geometry_model,
const GravityModel::Interface<dim> &gravity_model,
const MaterialModel::Interface<dim> &material_model,
+ const CompositionalInitialConditions::Interface<dim> &compositional_initial_conditions,
const double surface_pressure,
- const double surface_temperature)
+ const double surface_temperature,
+ const unsigned int n_compositional_fields)
:
n_points(1000),
temperatures(n_points, -1),
@@ -61,16 +63,22 @@
const Point<dim> representative_point = geometry_model.representative_point (z);
+ //TODO: we look up the composition at the representative point, but we should
+ // use averaged compositional values here. Right?
+ std::vector<double> initial_composition(n_compositional_fields);
+ for (unsigned int c=0;c<n_compositional_fields;++c)
+ initial_composition[c] = compositional_initial_conditions.initial_composition(representative_point, c);
+
// get material parameters and the magnitude of gravity. we assume
// that gravity always points along the depth direction. this
// may not strictly be true always but is likely a good enough
// approximation here.
const double density = material_model.density(temperatures[i-1], pressures[i-1],
- representative_point);
+ initial_composition,representative_point);
const double alpha = material_model.thermal_expansion_coefficient(temperatures[i-1], pressures[i-1],
- representative_point);
+ initial_composition,representative_point);
const double cp = material_model.specific_heat(temperatures[i-1], pressures[i-1],
- representative_point);
+ initial_composition, representative_point);
const double gravity = gravity_model.gravity_vector(representative_point).norm();
pressures[i] = pressures[i-1]
Property changes on: trunk/aspect/source/compositional_initial_conditions
___________________________________________________________________
Added: svn:ignore
+ Makefile.dep
Modified: trunk/aspect/source/compositional_initial_conditions/interface.cc
===================================================================
--- trunk/aspect/source/compositional_initial_conditions/interface.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/compositional_initial_conditions/interface.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -40,13 +40,9 @@
template <int dim>
void
- Interface<dim>::initialize (const GeometryModel::Interface<dim> &geometry_model_,
- const BoundaryTemperature::Interface<dim> &boundary_temperature_,
- const AdiabaticConditions<dim> &adiabatic_conditions_)
+ Interface<dim>::initialize (const GeometryModel::Interface<dim> &geometry_model_)
{
geometry_model = &geometry_model_;
- boundary_temperature = &boundary_temperature_;
- adiabatic_conditions = &adiabatic_conditions_;
}
@@ -94,9 +90,7 @@
template <int dim>
Interface<dim> *
create_initial_conditions (ParameterHandler &prm,
- const GeometryModel::Interface<dim> &geometry_model,
- const BoundaryTemperature::Interface<dim> &boundary_temperature,
- const AdiabaticConditions<dim> &adiabatic_conditions)
+ const GeometryModel::Interface<dim> &geometry_model)
{
// we need to get at the number of compositional fields here to
// know whether we need to generate something at all.
@@ -115,12 +109,10 @@
}
prm.leave_subsection ();
- Interface<dim> *plugin = std_cxx1x::get<dim>(registered_plugins).create_plugin (model_name, prm);
- plugin->initialize (geometry_model,
- boundary_temperature,
- adiabatic_conditions);
- return plugin;
- }
+ Interface<dim> *plugin = std_cxx1x::get<dim>(registered_plugins).create_plugin (model_name, prm);
+ plugin->initialize (geometry_model);
+ return plugin;
+ }
}
@@ -134,7 +126,7 @@
{
const std::string pattern_of_names
= std_cxx1x::get<dim>(registered_plugins).get_pattern_of_names ();
- prm.declare_entry ("Model name", "",
+ prm.declare_entry ("Model name", "function",
Patterns::Selection (pattern_of_names),
"Select one of the following models:\n\n"
+
@@ -182,9 +174,7 @@
template \
Interface<dim> * \
create_initial_conditions<dim> (ParameterHandler &prm, \
- const GeometryModel::Interface<dim> &geometry_model, \
- const BoundaryTemperature::Interface<dim> &boundary_temperature, \
- const AdiabaticConditions<dim> &adiabatic_conditions);
+ const GeometryModel::Interface<dim> &geometry_model);
ASPECT_INSTANTIATE(INSTANTIATE)
}
Modified: trunk/aspect/source/material_model/duretz_et_al.cc
===================================================================
--- trunk/aspect/source/material_model/duretz_et_al.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/material_model/duretz_et_al.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -38,6 +38,7 @@
SolCx<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &p) const
{
@@ -75,6 +76,7 @@
SolCx<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -93,6 +95,7 @@
SolCx<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -111,6 +114,7 @@
SolCx<dim>::
density (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &p) const
{
// defined as given in the paper, plus the constant
@@ -124,6 +128,7 @@
SolCx<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -135,6 +140,7 @@
SolCx<dim>::
compressibility (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0.0;
@@ -257,6 +263,7 @@
SolKz<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &p) const
{
@@ -295,6 +302,7 @@
SolKz<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -313,6 +321,7 @@
SolKz<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -331,6 +340,7 @@
SolKz<dim>::
density (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &p) const
{
// defined as given in the paper
@@ -343,6 +353,7 @@
SolKz<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -354,6 +365,7 @@
SolKz<dim>::
compressibility (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0.0;
@@ -420,6 +432,7 @@
Inclusion<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &p) const
{
@@ -457,6 +470,7 @@
Inclusion<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -475,6 +489,7 @@
Inclusion<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -493,6 +508,7 @@
Inclusion<dim>::
density (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &p) const
{
return 0;
@@ -504,6 +520,7 @@
Inclusion<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -515,6 +532,7 @@
Inclusion<dim>::
compressibility (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0.0;
Modified: trunk/aspect/source/material_model/interface.cc
===================================================================
--- trunk/aspect/source/material_model/interface.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/material_model/interface.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -21,6 +21,7 @@
#include <aspect/global.h>
+#include <aspect/simulator.h>
#include <aspect/material_model/interface.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/std_cxx1x/tuple.h>
@@ -47,6 +48,7 @@
double
Interface<dim>::viscosity_derivative (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &,
const NonlinearDependence::Dependence dependence) const
{
@@ -61,6 +63,7 @@
double
Interface<dim>::density_derivative (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &,
const NonlinearDependence::Dependence dependence) const
{
@@ -74,6 +77,7 @@
double
Interface<dim>::compressibility_derivative (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &,
const NonlinearDependence::Dependence dependence) const
{
@@ -87,6 +91,7 @@
double
Interface<dim>::specific_heat_derivative (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &,
const NonlinearDependence::Dependence dependence) const
{
@@ -100,6 +105,7 @@
double
Interface<dim>::thermal_conductivity_derivative (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &,
const NonlinearDependence::Dependence dependence) const
{
@@ -171,6 +177,7 @@
Interface<dim>::
viscosity_ratio (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const
{
@@ -183,6 +190,7 @@
Interface<dim>::
seismic_Vp (double dummy1,
double dummy2,
+ const std::vector<double> &, /*composition*/
const Point<dim> &dummy3) const
{
return -1.0;
@@ -194,6 +202,7 @@
Interface<dim>::
seismic_Vs (double dummy1,
double dummy2,
+ const std::vector<double> &, /*composition*/
const Point<dim> &dummy3) const
{
return -1.0;
@@ -204,7 +213,8 @@
unsigned int
Interface<dim>::
thermodynamic_phase (double dummy1,
- double dummy2) const
+ double dummy2,
+ const std::vector<double> & /*composition*/) const
{
return 0;
}
@@ -214,11 +224,12 @@
Interface<dim>::
thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
- return (-1./density(temperature, pressure, position)
+ return (-1./density(temperature, pressure, compositional_fields, position)
*
- density_derivative(temperature, pressure, position, NonlinearDependence::temperature));
+ density_derivative(temperature, pressure, compositional_fields, position, NonlinearDependence::temperature));
}
template <int dim>
@@ -241,6 +252,46 @@
std_cxx1x::get<dim>(registered_plugins).declare_parameters (prm);
}
+ template <int dim>
+ Interface<dim>::MaterialModelInputs::MaterialModelInputs(unsigned int n_points, unsigned int n_comp)
+ {
+ position.resize(n_points);
+ temperature.resize(n_points);
+ pressure.resize(n_points);
+ composition.resize(n_points);
+ for(unsigned int i=0;i<n_points;++i)
+ composition[i].resize(n_comp);
+ strain_rate.resize(n_points);
+ }
+
+ template <int dim>
+ Interface<dim>::MaterialModelOutputs::MaterialModelOutputs(unsigned int n_points)
+ {
+ viscosities.resize(n_points);
+ densities.resize(n_points);
+ thermal_expansion_coefficients.resize(n_points);
+ specific_heat.resize(n_points);
+ thermal_conductivities.resize(n_points);
+ compressibilities.resize(n_points);
+ }
+
+ template <int dim>
+ void
+ Interface<dim>::compute_parameters(const struct MaterialModelInputs &in, struct MaterialModelOutputs &out)
+ {
+ for (unsigned int i=0; i < in.temperature.size(); ++i)
+ {
+ out.viscosities[i] = viscosity (in.temperature[i], in.pressure[i], in.composition[i], in.strain_rate[i], in.position[i]);
+ out.densities[i] = density (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
+ out.thermal_expansion_coefficients[i] = thermal_expansion_coefficient (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
+ out.specific_heat[i] = specific_heat (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
+ out.thermal_conductivities[i] = thermal_conductivity (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
+ out.compressibilities[i] = compressibility (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
+ out.is_compressible = is_compressible();
+ }
+ }
+
+
}
}
Modified: trunk/aspect/source/material_model/simple.cc
===================================================================
--- trunk/aspect/source/material_model/simple.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/material_model/simple.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -35,6 +35,7 @@
Simple<dim>::
viscosity (const double temperature,
const double,
+ const std::vector<double> &composition, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &) const
{
@@ -42,6 +43,11 @@
const double temperature_dependence = std::max(std::min(std::exp(-thermal_viscosity_exponent*delta_temp/reference_T),1e2),1e-2);
return temperature_dependence * eta;
+/* return (this->n_compositional_fields()>0
+ ?
+ (6.5*composition[0]+1) * eta
+ :
+ eta);*/
}
@@ -74,6 +80,7 @@
Simple<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return reference_specific_heat;
@@ -92,6 +99,7 @@
Simple<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return k_value;
@@ -110,10 +118,15 @@
Simple<dim>::
density (const double temperature,
const double,
+ const std::vector<double> &compositional_fields, /*composition*/
const Point<dim> &) const
{
- return (reference_rho *
- (1 - thermal_alpha * (temperature - reference_T)));
+ return (this->n_compositional_fields()>0
+ ?
+ 100.0 * compositional_fields[0] + reference_rho *
+ (1 - thermal_alpha * (temperature - reference_T))
+ :
+ reference_rho * (1 - thermal_alpha * (temperature - reference_T)));
}
@@ -122,6 +135,7 @@
Simple<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return thermal_alpha;
@@ -133,13 +147,12 @@
Simple<dim>::
compressibility (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0.0;
}
-
-
template <int dim>
bool
Simple<dim>::
Modified: trunk/aspect/source/material_model/steinberger.cc
===================================================================
--- trunk/aspect/source/material_model/steinberger.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/material_model/steinberger.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -41,7 +41,8 @@
class material_lookup
{
public:
- material_lookup(const std::string &filename, const bool interpol)
+ material_lookup(const std::string &filename,
+ const bool interpol)
{
/* Initializing variables */
@@ -159,46 +160,65 @@
}
- double specific_heat(double temperature,double pressure) const
+ double
+ specific_heat(double temperature,
+ double pressure) const
{
return value(temperature,pressure,specific_heat_values,interpolation);
}
- double density(double temperature,double pressure) const
+ double
+ density(double temperature,
+ double pressure) const
{
return value(temperature,pressure,density_values,interpolation);
}
- double thermal_expansivity(const double temperature,const double pressure) const
+ double
+ thermal_expansivity(const double temperature,
+ const double pressure) const
{
return value(temperature,pressure,thermal_expansivity_values,interpolation);
}
- double seismic_Vp(const double temperature,const double pressure) const
+ double
+ seismic_Vp(const double temperature,
+ const double pressure) const
{
return value(temperature,pressure,vp_values,false);
}
- double seismic_Vs(const double temperature,const double pressure) const
+ double
+ seismic_Vs(const double temperature,
+ const double pressure) const
{
return value(temperature,pressure,vs_values,false);
}
- double dHdT (const double temperature,const double pressure) const
+ double
+ dHdT (const double temperature,
+ const double pressure) const
{
const double h = value(temperature,pressure,enthalpy_values,interpolation);
const double dh = value(temperature+delta_temp,pressure,enthalpy_values,interpolation);
return (dh - h) / delta_temp;
}
- double dHdp (const double temperature,const double pressure) const
+ double
+ dHdp (const double temperature,
+ const double pressure) const
{
const double h = value(temperature,pressure,enthalpy_values,interpolation);
const double dh = value(temperature,pressure+delta_press,enthalpy_values,interpolation);
return (dh - h) / delta_press;
}
- double value (const double temperature,const double pressure,const dealii::Table<2,double>& values,bool interpol) const
+ double
+ value (const double temperature,
+ const double pressure,
+ const dealii::Table<2,
+ double>& values,
+ bool interpol) const
{
const double nT = get_nT(temperature);
const unsigned int inT = static_cast<unsigned int>(nT);
@@ -391,6 +411,7 @@
Steinberger<dim>::
viscosity (const double temperature,
const double,
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &,
const Point<dim> &position) const
{
@@ -463,6 +484,7 @@
Steinberger<dim>::
specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
static std::string datapath = datadirectory+material_file_name;
@@ -483,6 +505,7 @@
Steinberger<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 4.7;
@@ -494,6 +517,7 @@
Steinberger<dim>::
density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
static std::string datapath = datadirectory+material_file_name;
@@ -505,6 +529,7 @@
Steinberger<dim>::
thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
static std::string datapath = datadirectory+material_file_name;
@@ -514,7 +539,7 @@
else
{
const double dHdp = get_material_data(datapath,interpolation).dHdp(temperature+get_deltat(position),pressure);
- const double alpha = (1 - density(temperature,pressure,position) * dHdp) / temperature;
+ const double alpha = (1 - density(temperature,pressure,compositional_fields,position) * dHdp) / temperature;
return std::max(std::min(alpha,1e-3),1e-5);
}
}
@@ -524,6 +549,7 @@
Steinberger<dim>::
seismic_Vp (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
static std::string datapath = datadirectory+material_file_name;
@@ -535,6 +561,7 @@
Steinberger<dim>::
seismic_Vs (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
static std::string datapath = datadirectory+material_file_name;
@@ -547,6 +574,7 @@
Steinberger<dim>::
compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const
{
return 0.0;
Modified: trunk/aspect/source/material_model/table.cc
===================================================================
--- trunk/aspect/source/material_model/table.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/material_model/table.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -347,6 +347,7 @@
Table<dim>::
viscosity (const double temperature,
const double pressure,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const
{
@@ -413,6 +414,7 @@
Table<dim>::
viscosity_ratio (const double temperature,
const double pressure,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &strain_rate,
const Point<dim> &position) const
{
@@ -466,6 +468,7 @@
Table<dim>::
thermal_expansion_coefficient (const double temperature,
const double pressure,
+ const std::vector<double> &, /*composition*/
const Point<dim> &p) const
{
static internal::P_T_LookupFunction alpha(data_directory+"alpha_bin");
@@ -477,6 +480,7 @@
Table<dim>::
specific_heat (const double temperature,
const double pressure,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
// const double reference_specific_heat = 1250; /* J / K / kg */ //??
@@ -500,6 +504,7 @@
Table<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
// this model assumes that the thermal conductivity is in fact constant
@@ -519,6 +524,7 @@
Table<dim>::
density (const double temperature,
const double pressure,
+ const std::vector<double> &, /*composition*/
const Point<dim> &position) const
{
static internal::P_T_LookupFunction rho(data_directory+"rho_bin");
@@ -529,7 +535,8 @@
double
Table<dim>::
seismic_Vp (const double temperature,
- const double pressure) const
+ const double pressure,
+ const std::vector<double> & /*composition*/) const
{
static internal::P_T_LookupFunction vp(data_directory+"vseis_p_bin");
return vp.value(temperature, pressure);
@@ -540,7 +547,8 @@
double
Table<dim>::
seismic_Vs (const double temperature,
- const double pressure) const
+ const double pressure,
+ const std::vector<double> & /*composition*/) const
{
static internal::P_T_LookupFunction vs(data_directory+"vseis_s_bin");
return vs.value(temperature, pressure);
@@ -551,7 +559,8 @@
unsigned int
Table<dim>::
thermodynamic_phase (const double temperature,
- const double pressure) const
+ const double pressure,
+ const std::vector<double> & /*composition*/) const
{
if (!compute_phases)
return 0;
@@ -566,6 +575,7 @@
Table<dim>::
compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &, /*composition*/
const Point<dim> &position) const
{
static internal::P_T_LookupFunction rho(data_directory+"rho_bin");
Modified: trunk/aspect/source/material_model/tan_gurnis.cc
===================================================================
--- trunk/aspect/source/material_model/tan_gurnis.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/material_model/tan_gurnis.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -50,6 +50,7 @@
TanGurnis<dim>::
viscosity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const SymmetricTensor<2,dim> &,
const Point<dim> &pos) const
{
@@ -87,6 +88,7 @@
TanGurnis<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 1250;
@@ -105,6 +107,7 @@
TanGurnis<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 2e-5;
@@ -123,6 +126,7 @@
TanGurnis<dim>::
density (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &pos) const
{
const double depth = 1.0-pos(dim-1);
@@ -136,6 +140,7 @@
TanGurnis<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return thermal_alpha;
@@ -147,9 +152,10 @@
TanGurnis<dim>::
compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &pos) const
{
- return Di/gamma / density(temperature, pressure, pos);
+ return Di/gamma / density(temperature, pressure, compositional_fields, pos);
}
Modified: trunk/aspect/source/postprocess/depth_average.cc
===================================================================
--- trunk/aspect/source/postprocess/depth_average.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/depth_average.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -75,8 +75,8 @@
// add temperature and the compositional fields that follow
// it immediately
this->get_depth_average_temperature(temp[i++]);
- for (unsigned int k=0; k<this->n_compositional_fields(); ++k)
- this->get_depth_average_composition(k, temp[i++]);
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+ this->get_depth_average_composition(c, temp[i++]);
this->get_adiabatic_conditions().get_adiabatic_temperature_profile(temp[i++],100);
this->get_depth_average_velocity_magnitude(temp[i++]);
this->get_depth_average_sinking_velocity(temp[i++]);
@@ -112,8 +112,8 @@
const std::string filename = this->get_output_directory() + "depthaverage.plt";
std::ofstream f (filename.c_str());
f << "# time, depth, avg T";
- for (unsigned int k=0; k<this->n_compositional_fields(); ++k)
- f << ", avg C_" << k+1;
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+ f << ", avg C_" << c+1;
f << ", adiabatic T, velocity magnitude, avg sinking velocity, avg Vs, avg Vp, avg viscosity" << std::endl;
f << std::scientific;
Modified: trunk/aspect/source/postprocess/heat_flux_statistics.cc
===================================================================
--- trunk/aspect/source/postprocess/heat_flux_statistics.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/heat_flux_statistics.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,10 +53,19 @@
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ compositional_fields.push_back(temp);
+ }
+
std::vector<Tensor<1,dim> > temperature_gradients (quadrature_formula.size());
std::vector<double> temperature_values (quadrature_formula.size());
std::vector<double> pressure_values (quadrature_formula.size());
+ std::vector<std::vector<double> > composition_values (this->n_compositional_fields(),std::vector<double> (quadrature_formula.size()));
+ std::vector<double> composition_values_at_q_point (this->n_compositional_fields());
std::map<types::boundary_id_t, double> local_boundary_fluxes;
@@ -85,13 +94,21 @@
temperature_values);
fe_face_values[pressure].get_function_values (this->get_solution(),
pressure_values);
+ for(unsigned int c=0;c<this->n_compositional_fields();++c)
+ fe_face_values[compositional_fields[c]].get_function_values(this->get_solution(),
+ composition_values[c]);
double local_normal_flux = 0;
for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
{
+ this->get_composition_values_at_q_point (composition_values,
+ q,
+ composition_values_at_q_point);
+
const double thermal_conductivity
= this->get_material_model().thermal_conductivity(temperature_values[q],
pressure_values[q],
+ composition_values_at_q_point,
fe_face_values.quadrature_point(q));
local_normal_flux
Modified: trunk/aspect/source/postprocess/table_heat_flux_statistics.cc
===================================================================
--- trunk/aspect/source/postprocess/table_heat_flux_statistics.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/table_heat_flux_statistics.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -52,11 +52,21 @@
update_q_points | update_JxW_values);
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ compositional_fields.push_back(temp);
+ }
+
std::vector<Tensor<1,dim> > temperature_gradients (quadrature_formula.size());
std::vector<double> temperature_values (quadrature_formula.size());
std::vector<double> pressure_values (quadrature_formula.size());
+ std::vector<std::vector<double> > composition_values (this->n_compositional_fields(),std::vector<double> (quadrature_formula.size()));
+ std::vector<double> composition_values_at_q_point (this->n_compositional_fields());
+
// find out which boundary indicators are related to Dirichlet temperature boundaries.
// it only makes sense to compute heat fluxes on these boundaries.
const std::set<types::boundary_id_t>
@@ -89,13 +99,21 @@
temperature_values);
fe_face_values[pressure].get_function_values (this->get_solution(),
pressure_values);
+ for(unsigned int c=0;c<this->n_compositional_fields();++c)
+ fe_face_values[compositional_fields[c]].get_function_values(this->get_solution(),
+ composition_values[c]);
double local_normal_flux = 0;
for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
{
+ this->get_composition_values_at_q_point (composition_values,
+ q,
+ composition_values_at_q_point);
+
const double thermal_conductivity
= this->get_material_model().thermal_conductivity(temperature_values[q],
pressure_values[q],
+ composition_values_at_q_point,
fe_face_values.quadrature_point(q));
local_normal_flux
Modified: trunk/aspect/source/postprocess/table_velocity_statistics.cc
===================================================================
--- trunk/aspect/source/postprocess/table_velocity_statistics.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/table_velocity_statistics.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -165,6 +165,8 @@
if (this->get_time() == 0e0)
{
double dT = this->get_boundary_temperature().maximal_temperature() - this->get_boundary_temperature().minimal_temperature();
+ // we do not compute the compositions but give the functions below the value 0.0 instead
+ std::vector<double> composition_values(this->n_compositional_fields(),0.0);
Point<dim> representative_point = Point<dim>::unit_vector(dim-1);
const double gravity = this->get_gravity_model().gravity_vector(representative_point).norm();
@@ -201,13 +203,13 @@
<< Ra
<< std::endl;
this->get_pcout()<< " k_value: "
- << material_model.thermal_conductivity(dT, dT, representative_point)
+ << material_model.thermal_conductivity(dT, dT, composition_values, representative_point) //TODO: dT for the pressure is wrong
<< std::endl;
this->get_pcout()<< " reference_cp: "
<< material_model.reference_cp()
<< std::endl;
this->get_pcout()<< " reference_thermal_diffusivity: "
- << material_model.thermal_conductivity(dT, dT, representative_point)/(material_model.reference_density()*material_model.reference_cp())
+ << material_model.thermal_conductivity(dT, dT, composition_values, representative_point)/(material_model.reference_density()*material_model.reference_cp()) //TODO: dT for the pressure is wrong
<< std::endl;
this->get_pcout()<< std::endl;
}
Modified: trunk/aspect/source/postprocess/velocity_statistics.cc
===================================================================
--- trunk/aspect/source/postprocess/velocity_statistics.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/velocity_statistics.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -148,6 +148,8 @@
const double h = this->get_geometry_model().maximal_depth();
const double dT = this->get_boundary_temperature().maximal_temperature() - this->get_boundary_temperature().minimal_temperature();
+ // we do not compute the compositions but give the functions below the value 0.0 instead
+ std::vector<double> composition_values(this->n_compositional_fields(),0.0);
Point<dim> representative_point = Point<dim>::unit_vector(dim-1);
const double gravity = this->get_gravity_model().gravity_vector(representative_point).norm();
@@ -184,13 +186,13 @@
<< Ra
<< std::endl;
this->get_pcout()<< " k_value: "
- << material_model.thermal_conductivity(dT, dT, representative_point)
+ << material_model.thermal_conductivity(dT, dT, composition_values, representative_point) //TODO: dT for the pressure is wrong
<< std::endl;
this->get_pcout()<< " reference_cp: "
<< material_model.reference_cp()
<< std::endl;
this->get_pcout()<< " reference_thermal_diffusivity: "
- << material_model.thermal_conductivity(dT, dT, representative_point)/(material_model.reference_density()*material_model.reference_cp())
+ << material_model.thermal_conductivity(dT, dT, composition_values, representative_point)/(material_model.reference_density()*material_model.reference_cp()) //TODO: dT for the pressure is wrong
<< std::endl;
this->get_pcout()<< std::endl;
}
Modified: trunk/aspect/source/postprocess/visualization/density.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/density.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/density.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,18 +53,22 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
// extract the primal variables
const double pressure = uh[q][dim];
const double temperature = uh[q][dim+1];
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ composition[c] = uh[q][dim+2+c];
computed_quantities[q](0) = this->get_material_model().density(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: trunk/aspect/source/postprocess/visualization/friction_heating.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/friction_heating.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/friction_heating.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,10 +53,10 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
- Assert (duh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (duh[0].size() == dim+2+this->n_compositional_fields(),ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
@@ -75,8 +75,14 @@
strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
:
strain_rate);
+
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ composition[c] = uh[q][dim+2+c];
+
computed_quantities[q](0) = 2 * this->get_material_model().viscosity(temperature,
pressure,
+ composition,
strain_rate,
evaluation_points[q]) *
compressible_strain_rate * compressible_strain_rate;
Modified: trunk/aspect/source/postprocess/visualization/nonadiabatic_pressure.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/nonadiabatic_pressure.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/nonadiabatic_pressure.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,9 +53,9 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
Modified: trunk/aspect/source/postprocess/visualization/nonadiabatic_temperature.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/nonadiabatic_temperature.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/nonadiabatic_temperature.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,9 +53,9 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
Modified: trunk/aspect/source/postprocess/visualization/seismic_vp.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/seismic_vp.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/seismic_vp.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,17 +53,21 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
const double pressure = uh[q][dim];
const double temperature = uh[q][dim+1];
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ composition[c] = uh[q][dim+2+c];
computed_quantities[q](0) = this->get_material_model().seismic_Vp(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: trunk/aspect/source/postprocess/visualization/seismic_vs.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/seismic_vs.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/seismic_vs.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,17 +53,21 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
const double pressure = uh[q][dim];
const double temperature = uh[q][dim+1];
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ composition[c] = uh[q][dim+2+c];
computed_quantities[q](0) = this->get_material_model().seismic_Vs(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: trunk/aspect/source/postprocess/visualization/specific_heat.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/specific_heat.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/specific_heat.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,17 +53,21 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
const double pressure = uh[q][dim];
const double temperature = uh[q][dim+1];
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ composition[c] = uh[q][dim+2+c];
computed_quantities[q](0) = this->get_material_model().specific_heat(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: trunk/aspect/source/postprocess/visualization/strain_rate.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/strain_rate.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/strain_rate.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,10 +53,10 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
- Assert (duh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (duh[0].size() == dim+2+this->n_compositional_fields(),ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
Modified: trunk/aspect/source/postprocess/visualization/thermal_expansivity.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/thermal_expansivity.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/thermal_expansivity.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,18 +53,22 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
// extract the primal variables
const double pressure = uh[q][dim];
const double temperature = uh[q][dim+1];
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ composition[c] = uh[q][dim+2+c];
computed_quantities[q](0) = this->get_material_model().thermal_expansion_coefficient(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: trunk/aspect/source/postprocess/visualization/thermodynamic_phase.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/thermodynamic_phase.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/thermodynamic_phase.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,17 +53,21 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
const double pressure = uh[q][dim];
const double temperature = uh[q][dim+1];
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ composition[c] = uh[q][dim+2+c];
computed_quantities[q](0) = this->get_material_model().thermodynamic_phase(temperature,
- pressure);
+ pressure,
+ composition);
}
}
}
Modified: trunk/aspect/source/postprocess/visualization/viscosity.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/viscosity.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/viscosity.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,10 +53,10 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
- Assert (duh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (duh[0].size() == dim+2+this->n_compositional_fields(),ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
@@ -75,8 +75,14 @@
strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
:
strain_rate);
+
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int i=0;i<this->n_compositional_fields();++i)
+ composition[i] = uh[q][dim+2+i];
+
computed_quantities[q](0) = this->get_material_model().viscosity(temperature,
pressure,
+ composition,
strain_rate,
evaluation_points[q]);
}
Modified: trunk/aspect/source/postprocess/visualization/viscosity_ratio.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization/viscosity_ratio.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization/viscosity_ratio.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -53,10 +53,10 @@
std::vector<Vector<double> > &computed_quantities) const
{
const unsigned int n_quadrature_points = uh.size();
- Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
- Assert (computed_quantities[0].size() == 1, ExcInternalError());
- Assert (uh[0].size() == dim+2, ExcInternalError());
- Assert (duh[0].size() == dim+2, ExcInternalError());
+ Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError());
+ Assert (computed_quantities[0].size() == 1, ExcInternalError());
+ Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+ Assert (duh[0].size() == dim+2+this->n_compositional_fields(),ExcInternalError());
for (unsigned int q=0; q<n_quadrature_points; ++q)
{
@@ -75,8 +75,14 @@
strain_rate - 1./3 * trace(strain_rate) * unit_symmetric_tensor<dim>()
:
strain_rate);
+
+ std::vector<double> composition(this->n_compositional_fields());
+ for (unsigned int c=0;c<this->n_compositional_fields();++c)
+ composition[c] = uh[q][dim+2+c];
+
computed_quantities[q](0) = this->get_material_model().viscosity_ratio(temperature,
pressure,
+ composition,
strain_rate,
evaluation_points[q]);
}
Modified: trunk/aspect/source/postprocess/visualization.cc
===================================================================
--- trunk/aspect/source/postprocess/visualization.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/postprocess/visualization.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -129,8 +129,8 @@
std::vector<std::string> solution_names (dim, "velocity");
solution_names.push_back ("p");
solution_names.push_back ("T");
- for (unsigned int i=0; i<this->n_compositional_fields(); ++i)
- solution_names.push_back ("C_" + boost::lexical_cast<std::string>(i+1));
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+ solution_names.push_back ("C_" + boost::lexical_cast<std::string>(c+1));
std::vector<DataComponentInterpretation::DataComponentInterpretation>
@@ -138,7 +138,7 @@
DataComponentInterpretation::component_is_part_of_vector);
interpretation.push_back (DataComponentInterpretation::component_is_scalar);
interpretation.push_back (DataComponentInterpretation::component_is_scalar);
- for (unsigned int i=0; i<this->n_compositional_fields(); ++i)
+ for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
interpretation.push_back (DataComponentInterpretation::component_is_scalar);
data_out.add_data_vector (this->get_solution(),
Modified: trunk/aspect/source/simulator/assembly.cc
===================================================================
--- trunk/aspect/source/simulator/assembly.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/simulator/assembly.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -49,7 +49,8 @@
StokesPreconditioner (const FiniteElement<dim> &finite_element,
const Quadrature<dim> &quadrature,
const Mapping<dim> &mapping,
- const UpdateFlags update_flags);
+ const UpdateFlags update_flags,
+ const unsigned int n_compositional_fields);
StokesPreconditioner (const StokesPreconditioner &data);
FEValues<dim> finite_element_values;
@@ -60,6 +61,10 @@
std::vector<double> temperature_values;
std::vector<double> pressure_values;
std::vector<SymmetricTensor<2,dim> > strain_rates;
+ std::vector<std::vector<double> > composition_values;
+
+ typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs material_model_outputs;
};
@@ -69,7 +74,8 @@
StokesPreconditioner (const FiniteElement<dim> &finite_element,
const Quadrature<dim> &quadrature,
const Mapping<dim> &mapping,
- const UpdateFlags update_flags)
+ const UpdateFlags update_flags,
+ const unsigned int n_compositional_fields)
:
finite_element_values (mapping, finite_element, quadrature,
update_flags),
@@ -77,7 +83,11 @@
phi_p (finite_element.dofs_per_cell),
temperature_values (quadrature.size()),
pressure_values (quadrature.size()),
- strain_rates (quadrature.size())
+ strain_rates (quadrature.size()),
+ composition_values(n_compositional_fields,
+ std::vector<double>(quadrature.size())),
+ material_model_inputs(quadrature.size(), n_compositional_fields),
+ material_model_outputs(quadrature.size())
{}
@@ -94,7 +104,10 @@
phi_p (scratch.phi_p),
temperature_values (scratch.temperature_values),
pressure_values (scratch.pressure_values),
- strain_rates (scratch.strain_rates)
+ strain_rates (scratch.strain_rates),
+ composition_values(scratch.composition_values),
+ material_model_inputs(scratch.material_model_inputs),
+ material_model_outputs(scratch.material_model_outputs)
{}
@@ -120,7 +133,8 @@
StokesSystem (const FiniteElement<dim> &finite_element,
const Mapping<dim> &mapping,
const Quadrature<dim> &quadrature,
- const UpdateFlags update_flags);
+ const UpdateFlags update_flags,
+ const unsigned int n_compositional_fields);
StokesSystem (const StokesSystem<dim> &data);
@@ -128,6 +142,10 @@
std::vector<SymmetricTensor<2,dim> > grads_phi_u;
std::vector<double> div_phi_u;
std::vector<Tensor<1,dim> > velocity_values;
+ std::vector<std::vector<double> > composition_values;
+
+ typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs material_model_outputs;
};
@@ -137,15 +155,20 @@
StokesSystem (const FiniteElement<dim> &finite_element,
const Mapping<dim> &mapping,
const Quadrature<dim> &quadrature,
- const UpdateFlags update_flags)
+ const UpdateFlags update_flags,
+ const unsigned int n_compositional_fields)
:
StokesPreconditioner<dim> (finite_element, quadrature,
mapping,
- update_flags),
+ update_flags, n_compositional_fields),
phi_u (finite_element.dofs_per_cell),
grads_phi_u (finite_element.dofs_per_cell),
div_phi_u (finite_element.dofs_per_cell),
- velocity_values (quadrature.size())
+ velocity_values (quadrature.size()),
+ composition_values(n_compositional_fields,
+ std::vector<double>(quadrature.size())),
+ material_model_inputs(quadrature.size(), n_compositional_fields),
+ material_model_outputs(quadrature.size())
{}
@@ -158,7 +181,10 @@
phi_u (scratch.phi_u),
grads_phi_u (scratch.grads_phi_u),
div_phi_u (scratch.div_phi_u),
- velocity_values (scratch.velocity_values)
+ velocity_values (scratch.velocity_values),
+ composition_values(scratch.composition_values),
+ material_model_inputs(scratch.material_model_inputs),
+ material_model_outputs(scratch.material_model_outputs)
{}
@@ -168,7 +194,8 @@
{
TemperatureSystem (const FiniteElement<dim> &finite_element,
const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature);
+ const Quadrature<dim> &quadrature,
+ const unsigned int n_compositional_fields);
TemperatureSystem (const TemperatureSystem &data);
FEValues<dim> finite_element_values;
@@ -192,10 +219,18 @@
std::vector<double> old_temperature_laplacians;
std::vector<double> old_old_temperature_laplacians;
+ std::vector<std::vector<double> > old_composition_values;
+ std::vector<std::vector<double> > old_old_composition_values;
+
std::vector<double> current_temperature_values;
std::vector<Tensor<1,dim> > current_velocity_values;
std::vector<SymmetricTensor<2,dim> > current_strain_rates;
std::vector<double> current_pressure_values;
+ std::vector<std::vector<double> > current_composition_values;
+
+
+ typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs;
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs material_model_outputs;
};
@@ -204,7 +239,8 @@
TemperatureSystem<dim>::
TemperatureSystem (const FiniteElement<dim> &finite_element,
const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature)
+ const Quadrature<dim> &quadrature,
+ const unsigned int n_compositional_fields)
:
finite_element_values (mapping,
finite_element, quadrature,
@@ -228,10 +264,18 @@
old_old_temperature_grads(quadrature.size()),
old_temperature_laplacians(quadrature.size()),
old_old_temperature_laplacians(quadrature.size()),
+ old_composition_values(n_compositional_fields,
+ std::vector<double>(quadrature.size())),
+ old_old_composition_values(n_compositional_fields,
+ std::vector<double>(quadrature.size())),
current_temperature_values(quadrature.size()),
current_velocity_values(quadrature.size()),
current_strain_rates(quadrature.size()),
- current_pressure_values(quadrature.size())
+ current_pressure_values(quadrature.size()),
+ current_composition_values(n_compositional_fields,
+ std::vector<double>(quadrature.size())),
+ material_model_inputs(quadrature.size(), n_compositional_fields),
+ material_model_outputs(quadrature.size())
{}
@@ -259,10 +303,15 @@
old_old_temperature_grads (scratch.old_old_temperature_grads),
old_temperature_laplacians (scratch.old_temperature_laplacians),
old_old_temperature_laplacians (scratch.old_old_temperature_laplacians),
+ old_composition_values(scratch.old_composition_values),
+ old_old_composition_values(scratch.old_old_composition_values),
current_temperature_values(scratch.current_temperature_values),
current_velocity_values(scratch.current_velocity_values),
current_strain_rates(scratch.current_strain_rates),
- current_pressure_values(scratch.current_pressure_values)
+ current_pressure_values(scratch.current_pressure_values),
+ current_composition_values(scratch.current_composition_values),
+ material_model_inputs(scratch.material_model_inputs),
+ material_model_outputs(scratch.material_model_outputs)
{}
@@ -271,7 +320,8 @@
{
CompositionSystem (const FiniteElement<dim> &finite_element,
const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature);
+ const Quadrature<dim> &quadrature,
+ const unsigned int n_compositional_fields);
CompositionSystem (const CompositionSystem &data);
FEValues<dim> finite_element_values;
@@ -282,12 +332,6 @@
std::vector<Tensor<1,dim> > old_velocity_values;
std::vector<Tensor<1,dim> > old_old_velocity_values;
- std::vector<double> old_pressure;
- std::vector<double> old_old_pressure;
-
- std::vector<SymmetricTensor<2,dim> > old_strain_rates;
- std::vector<SymmetricTensor<2,dim> > old_old_strain_rates;
-
std::vector<double> old_composition_values;
std::vector<double> old_old_composition_values;
std::vector<Tensor<1,dim> > old_composition_grads;
@@ -295,10 +339,7 @@
std::vector<double> old_composition_laplacians;
std::vector<double> old_old_composition_laplacians;
- std::vector<double> current_composition_values;
std::vector<Tensor<1,dim> > current_velocity_values;
- std::vector<SymmetricTensor<2,dim> > current_strain_rates;
- std::vector<double> current_pressure_values;
};
@@ -307,7 +348,8 @@
CompositionSystem<dim>::
CompositionSystem (const FiniteElement<dim> &finite_element,
const Mapping<dim> &mapping,
- const Quadrature<dim> &quadrature)
+ const Quadrature<dim> &quadrature,
+ const unsigned int n_compositional_fields)
:
finite_element_values (mapping,
finite_element, quadrature,
@@ -321,20 +363,13 @@
grad_phi_C (finite_element.dofs_per_cell),
old_velocity_values (quadrature.size()),
old_old_velocity_values (quadrature.size()),
- old_pressure (quadrature.size()),
- old_old_pressure (quadrature.size()),
- old_strain_rates (quadrature.size()),
- old_old_strain_rates (quadrature.size()),
- old_composition_values (quadrature.size()),
+ old_composition_values(quadrature.size()),
old_old_composition_values(quadrature.size()),
old_composition_grads(quadrature.size()),
old_old_composition_grads(quadrature.size()),
old_composition_laplacians(quadrature.size()),
old_old_composition_laplacians(quadrature.size()),
- current_composition_values(quadrature.size()),
- current_velocity_values(quadrature.size()),
- current_strain_rates(quadrature.size()),
- current_pressure_values(quadrature.size())
+ current_velocity_values(quadrature.size())
{}
@@ -352,20 +387,13 @@
grad_phi_C (scratch.grad_phi_C),
old_velocity_values (scratch.old_velocity_values),
old_old_velocity_values (scratch.old_old_velocity_values),
- old_pressure (scratch.old_pressure),
- old_old_pressure (scratch.old_old_pressure),
- old_strain_rates (scratch.old_strain_rates),
- old_old_strain_rates (scratch.old_old_strain_rates),
old_composition_values (scratch.old_composition_values),
old_old_composition_values (scratch.old_old_composition_values),
old_composition_grads (scratch.old_composition_grads),
old_old_composition_grads (scratch.old_old_composition_grads),
old_composition_laplacians (scratch.old_composition_laplacians),
old_old_composition_laplacians (scratch.old_old_composition_laplacians),
- current_composition_values(scratch.current_composition_values),
- current_velocity_values(scratch.current_velocity_values),
- current_strain_rates(scratch.current_strain_rates),
- current_pressure_values(scratch.current_pressure_values)
+ current_velocity_values(scratch.current_velocity_values)
{}
}
@@ -624,70 +652,72 @@
}
-
template <int dim>
- double
+ void
Simulator<dim>::
- compute_viscosity (const std::vector<double> &old_temperature,
- const std::vector<double> &old_old_temperature,
- const std::vector<Tensor<1,dim> > &old_temperature_grads,
- const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
- const std::vector<double> &old_temperature_laplacians,
- const std::vector<double> &old_old_temperature_laplacians,
- const std::vector<Tensor<1,dim> > &old_velocity_values,
- const std::vector<Tensor<1,dim> > &old_old_velocity_values,
- const std::vector<SymmetricTensor<2,dim> > &old_strain_rates,
- const std::vector<SymmetricTensor<2,dim> > &old_old_strain_rates,
- const std::vector<double> &old_pressure,
- const std::vector<double> &old_old_pressure,
- const double global_u_infty,
- const double global_T_variation,
- const double average_temperature,
- const double global_entropy_variation,
- const std::vector<Point<dim> > &evaluation_points,
- const double cell_diameter) const
- {
- if (global_u_infty == 0)
- return 5e-3 * cell_diameter;
+ compute_temperature_system_residual(const std::vector<double> &old_temperature,
+ const std::vector<double> &old_old_temperature,
+ const std::vector<Tensor<1,dim> > &old_temperature_grads,
+ const std::vector<Tensor<1,dim> > &old_old_temperature_grads,
+ const std::vector<double> &old_temperature_laplacians,
+ const std::vector<double> &old_old_temperature_laplacians,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
+ const std::vector<SymmetricTensor<2,dim> > &old_strain_rates,
+ const std::vector<SymmetricTensor<2,dim> > &old_old_strain_rates,
+ const std::vector<double> &old_pressure,
+ const std::vector<double> &old_old_pressure,
+ const std::vector<std::vector<double> > &old_composition,
+ const std::vector<std::vector<double> > &old_old_composition,
+ const double average_temperature,
+ const std::vector<Point<dim> > &evaluation_points,
+ double &max_residual,
+ double &max_velocity) const
+ {
const unsigned int n_q_points = old_temperature.size();
- double max_residual = 0;
- double max_velocity = 0;
+ typename MaterialModel::Interface<dim>::MaterialModelInputs inputs (n_q_points, parameters.n_compositional_fields);
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs outputs (n_q_points);
+ for (unsigned int q=0; q<n_q_points; ++q) {
+ inputs.temperature[q] = (old_temperature[q] + old_old_temperature[q]) / 2;
+ inputs.position[q] = evaluation_points[q];
+ inputs.pressure[q] = (old_pressure[q] + old_old_pressure[q]) / 2;
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+ inputs.composition[q][c] = (old_composition[c][q] + old_old_composition[c][q]) / 2;
+ inputs.strain_rate[q] = (old_strain_rates[q] + old_old_strain_rates[q]) / 2;
+ }
+
+ material_model->compute_parameters(inputs,outputs);
+
for (unsigned int q=0; q < n_q_points; ++q)
{
const Tensor<1,dim> u = (old_velocity_values[q] +
old_old_velocity_values[q]) / 2;
- const SymmetricTensor<2,dim> strain_rate = (old_strain_rates[q] +
- old_old_strain_rates[q]) / 2;
+ const SymmetricTensor<2,dim> strain_rate = inputs.strain_rate[q];
- const double T = (old_temperature[q] + old_old_temperature[q]) / 2;
- const double p = (old_pressure[q] + old_old_pressure[q]) / 2;
const double dT_dt = (old_temperature[q] - old_old_temperature[q])
/ old_time_step;
const double u_grad_T = u * (old_temperature_grads[q] +
old_old_temperature_grads[q]) / 2;
- const double alpha = material_model->thermal_expansion_coefficient(T, p, evaluation_points[q]);
- const double density = material_model->density(T, p, evaluation_points[q]);
- const double thermal_conductivity = material_model->thermal_conductivity(T, p, evaluation_points[q]);
- const double c_P = material_model->specific_heat(T, p, evaluation_points[q]);
+ const double alpha = outputs.thermal_expansion_coefficients[q];
+ const double density = outputs.densities[q];
+ const double thermal_conductivity = outputs.thermal_conductivities[q];
+ const double c_P = outputs.specific_heat[q];
const double k_Delta_T = thermal_conductivity
* (old_temperature_laplacians[q] +
old_old_temperature_laplacians[q]) / 2;
// verify correctness of the heating term
- const double viscosity = material_model->viscosity(T,
- p,
- strain_rate,
- evaluation_points[q]);
- const bool is_compressible = material_model->is_compressible ();
+ const double viscosity = outputs.viscosities[q];
+ const bool is_compressible = outputs.is_compressible;
const double compressibility
= (is_compressible
?
- material_model->compressibility(T, p, evaluation_points[q] )
+ outputs.compressibilities[q]
:
std::numeric_limits<double>::quiet_NaN() );
const Tensor<1,dim> gravity = gravity_model->gravity_vector (evaluation_points[q] );
@@ -728,19 +758,143 @@
// alpha = - 1/rho drho/dT
(parameters.include_adiabatic_heating
?
- alpha * density * (u*gravity) * T
+ alpha * density * (u*gravity) * inputs.temperature[q]
:
0)
);
double residual
= std::abs(density * c_P * (dT_dt + u_grad_T) - k_Delta_T - gamma);
if (parameters.stabilization_alpha == 2)
- residual *= std::abs(T - average_temperature);
+ residual *= std::abs(inputs.temperature[q] - average_temperature);
max_residual = std::max (residual, max_residual);
max_velocity = std::max (std::sqrt (u*u), max_velocity);
}
+ }
+ template <int dim>
+ void
+ Simulator<dim>::
+ compute_composition_system_residual(const std::vector<double> &old_composition,
+ const std::vector<double> &old_old_composition,
+ const std::vector<Tensor<1,dim> > &old_composition_grads,
+ const std::vector<Tensor<1,dim> > &old_old_composition_grads,
+ const std::vector<double> &old_composition_laplacians,
+ const std::vector<double> &old_old_composition_laplacians,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
+ const double average_composition,
+ const unsigned int composition_index,
+ double &max_residual,
+ double &max_velocity) const
+ {
+
+ const unsigned int n_q_points = old_composition.size();
+
+ for (unsigned int q=0; q < n_q_points; ++q)
+ {
+ const Tensor<1,dim> u = (old_velocity_values[q] +
+ old_old_velocity_values[q]) / 2;
+ const double C = (old_composition[q] + old_old_composition[q]) / 2;
+
+ const double dC_dt = (old_composition[q] - old_old_composition[q])
+ / old_time_step;
+ const double u_grad_C = u * (old_composition_grads[q] +
+ old_old_composition_grads[q]) / 2;
+
+ // the chemical diffusivity kappa is always smaller than the numerical diffusivity
+ // therefore we set it to zero
+ const double kappa = 0;
+
+ const double kappa_Delta_C = kappa
+ * (old_composition_laplacians[q] +
+ old_old_composition_laplacians[q]) / 2;
+
+ double residual
+ = std::abs(dC_dt + u_grad_C - kappa_Delta_C);
+ if (parameters.stabilization_alpha == 2)
+ residual *= std::abs(C - average_composition);
+
+ max_residual = std::max (residual, max_residual);
+ max_velocity = std::max (std::sqrt (u*u), max_velocity);
+ }
+ }
+
+
+
+ template <int dim>
+ double
+ Simulator<dim>::
+ compute_viscosity (const std::vector<double> &old_field,
+ const std::vector<double> &old_old_field,
+ const std::vector<Tensor<1,dim> > &old_field_grads,
+ const std::vector<Tensor<1,dim> > &old_old_field_grads,
+ const std::vector<double> &old_field_laplacians,
+ const std::vector<double> &old_old_field_laplacians,
+ const std::vector<Tensor<1,dim> > &old_velocity_values,
+ const std::vector<Tensor<1,dim> > &old_old_velocity_values,
+ const std::vector<SymmetricTensor<2,dim> > &old_strain_rates,
+ const std::vector<SymmetricTensor<2,dim> > &old_old_strain_rates,
+ const std::vector<double> &old_pressure,
+ const std::vector<double> &old_old_pressure,
+ const std::vector<std::vector<double> > &old_composition,
+ const std::vector<std::vector<double> > &old_old_composition,
+ const double global_u_infty,
+ const double global_T_variation,
+ const double average_field,
+ const double global_entropy_variation,
+ const std::vector<Point<dim> > &evaluation_points,
+ const double cell_diameter,
+ const unsigned int index) const
+ {
+ if (global_u_infty == 0)
+ return 5e-3 * cell_diameter;
+
+ double max_residual = 0;
+ double max_velocity = 0;
+
+ AssertIndexRange(index,parameters.n_compositional_fields+1);
+
+ if(index == 0) {
+ //make sure that all arguments we need for computing the residual are passed
+ Assert (old_strain_rates.size() > 0 && old_old_strain_rates.size() > 0
+ && old_pressure.size() > 0 && old_old_pressure.size() > 0,
+ ExcMessage ("Not enough parameters to calculate artificial viscosity "
+ "for the temperature equation."));
+
+ compute_temperature_system_residual(old_field,
+ old_old_field,
+ old_field_grads,
+ old_old_field_grads,
+ old_field_laplacians,
+ old_old_field_laplacians,
+ old_velocity_values,
+ old_old_velocity_values,
+ old_strain_rates,
+ old_old_strain_rates,
+ old_pressure,
+ old_old_pressure,
+ old_composition,
+ old_old_composition,
+ average_field,
+ evaluation_points,
+ max_residual,
+ max_velocity);
+ }
+ else
+ compute_composition_system_residual(old_field,
+ old_old_field,
+ old_field_grads,
+ old_old_field_grads,
+ old_field_laplacians,
+ old_old_field_laplacians,
+ old_velocity_values,
+ old_old_velocity_values,
+ average_field,
+ index-1, //index of compositional field
+ max_residual,
+ max_velocity);
+
const double max_viscosity = (parameters.stabilization_beta *
max_velocity * cell_diameter);
if (timestep_number <= 1)
@@ -767,7 +921,54 @@
}
+ template <int dim>
+ void
+ Simulator<dim>::
+ compute_material_model_input_values (const TrilinosWrappers::MPI::BlockVector ¤t_linearization_point,
+ const FEValues<dim> &input_finite_element_values,
+ const bool compute_strainrate,
+ typename MaterialModel::Interface<dim>::MaterialModelInputs &material_model_inputs) const
+ {
+ const FEValuesExtractors::Vector input_velocities (0);
+ const FEValuesExtractors::Scalar input_pressure (dim);
+ const FEValuesExtractors::Scalar input_temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> input_composition;
+
+ unsigned int n_q_points = material_model_inputs.temperature.size();
+
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ input_composition.push_back(temp);
+ }
+
+ input_finite_element_values[input_temperature].get_function_values (current_linearization_point,
+ material_model_inputs.temperature);
+ input_finite_element_values[input_pressure].get_function_values(current_linearization_point,
+ material_model_inputs.pressure);
+ if(compute_strainrate)
+ input_finite_element_values[input_velocities].get_function_symmetric_gradients(current_linearization_point,
+ material_model_inputs.strain_rate);
+
+ // the values of the compositional fields are stored as blockvectors for each field
+ // we have to extract them in this structure
+ std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,
+ std::vector<double> (n_q_points));
+
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ input_finite_element_values[input_composition[c]].get_function_values(current_linearization_point,
+ composition_values[c]);
+
+ // then we copy these values to exchange the inner and outer vector, because for the material
+ // model we need a vector with values of all the compositional fields for every quadrature point
+ for(unsigned int q=0;q<n_q_points;++q)
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ material_model_inputs.composition[q][c] = composition_values[c][q];
+
+ }
+
+
template <int dim>
void
Simulator<dim>::
@@ -780,23 +981,22 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
- const FEValuesExtractors::Scalar temperature (dim+1);
scratch.finite_element_values.reinit (cell);
- scratch.finite_element_values[temperature].get_function_values (current_linearization_point,
- scratch.temperature_values);
- scratch.finite_element_values[pressure].get_function_values(current_linearization_point,
- scratch.pressure_values);
- scratch.finite_element_values[velocities].get_function_symmetric_gradients(current_linearization_point,
- scratch.strain_rates);
-
data.local_matrix = 0;
+ compute_material_model_input_values (current_linearization_point,
+ scratch.finite_element_values,
+ true,
+ scratch.material_model_inputs);
+
+ material_model->compute_parameters(scratch.material_model_inputs,scratch.material_model_outputs);
+
for (unsigned int q=0; q<n_q_points; ++q)
{
- const double current_temperature = scratch.temperature_values[q];
- const double current_pressure = scratch.pressure_values[q];
+ const double current_temperature = scratch.material_model_inputs.temperature[q];
+ const double current_pressure = scratch.material_model_inputs.pressure[q];
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
@@ -804,10 +1004,7 @@
scratch.phi_p[k] = scratch.finite_element_values[pressure].value (k, q);
}
- const double eta = material_model->viscosity(current_temperature,
- current_pressure,
- scratch.strain_rates[q],
- scratch.finite_element_values.quadrature_point(q) );
+ const double eta = scratch.material_model_outputs.viscosities[q];
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
@@ -875,7 +1072,8 @@
update_JxW_values |
update_values |
update_gradients |
- update_quadrature_points),
+ update_quadrature_points,
+ parameters.n_compositional_fields),
internal::Assembly::CopyData::
StokesPreconditioner<dim> (finite_element));
@@ -940,42 +1138,31 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
- const FEValuesExtractors::Scalar temperature (dim+1);
scratch.finite_element_values.reinit (cell);
- //scratch.finite_element_values[temperature].get_function_values (old_solution,
- // scratch.old_temperature_values);
- // Assuming we already have the temperature for the current time step:
- scratch.finite_element_values[temperature].get_function_values (current_linearization_point,
- scratch.temperature_values);
- scratch.finite_element_values[pressure].get_function_values(current_linearization_point,
- scratch.pressure_values);
- scratch.finite_element_values[velocities].get_function_values(current_linearization_point,
- scratch.velocity_values);
- // we only need the strain rates for the viscosity,
- // which we only need when rebuilding the matrix
if (rebuild_stokes_matrix)
- scratch.finite_element_values[velocities]
- .get_function_symmetric_gradients(current_linearization_point,
- scratch.strain_rates);
-
-
-
- // cache whether the model is compressible or not
- const bool is_compressible = material_model->is_compressible ();
-
- if (rebuild_stokes_matrix)
data.local_matrix = 0;
data.local_rhs = 0;
- if (is_compressible)
+ if (material_model->is_compressible())
data.local_pressure_shape_function_integrals = 0;
+ // we only need the strain rates for the viscosity,
+ // which we only need when rebuilding the matrix
+ compute_material_model_input_values (current_linearization_point,
+ scratch.finite_element_values,
+ rebuild_stokes_matrix,
+ scratch.material_model_inputs);
+ material_model->compute_parameters(scratch.material_model_inputs,scratch.material_model_outputs);
+
+ scratch.finite_element_values[velocities].get_function_values(current_linearization_point,
+ scratch.velocity_values);
+
for (unsigned int q=0; q<n_q_points; ++q)
{
- const double current_temperature = scratch.temperature_values[q];
- const double current_pressure = scratch.pressure_values[q];
+ const double current_temperature = scratch.material_model_inputs.temperature[q];
+ const double current_pressure = scratch.material_model_inputs.pressure[q];
const Tensor<1,dim> current_u = scratch.velocity_values[q];
for (unsigned int k=0; k<dofs_per_cell; ++k)
@@ -991,10 +1178,7 @@
const double eta = (rebuild_stokes_matrix
?
- material_model->viscosity(current_temperature,
- current_pressure,
- scratch.strain_rates[q],
- scratch.finite_element_values.quadrature_point(q))
+ scratch.material_model_outputs.viscosities[q]
:
std::numeric_limits<double>::quiet_NaN());
@@ -1002,24 +1186,18 @@
gravity = gravity_model->gravity_vector (scratch.finite_element_values.quadrature_point(q));
const double compressibility
- = (is_compressible
+ = (scratch.material_model_outputs.is_compressible
?
- material_model->compressibility(current_temperature,
- current_pressure,
- scratch.finite_element_values
- .quadrature_point(q))
+ scratch.material_model_outputs.compressibilities[q]
:
std::numeric_limits<double>::quiet_NaN() );
- const double density = material_model->density(current_temperature,
- current_pressure,
- scratch.finite_element_values
- .quadrature_point(q));
+ const double density = scratch.material_model_outputs.densities[q];
if (rebuild_stokes_matrix)
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
data.local_matrix(i,j) += ( eta * 2.0 * (scratch.grads_phi_u[i] * scratch.grads_phi_u[j])
- - (is_compressible
+ - (scratch.material_model_outputs.is_compressible
?
eta * 2.0/3.0 * (scratch.div_phi_u[i] * scratch.div_phi_u[j])
:
@@ -1033,7 +1211,7 @@
for (unsigned int i=0; i<dofs_per_cell; ++i)
data.local_rhs(i) += (
(density * gravity * scratch.phi_u[i])
- + (is_compressible
+ + (scratch.material_model_outputs.is_compressible
?
(pressure_scaling *
compressibility * density *
@@ -1043,7 +1221,7 @@
0)
)
* scratch.finite_element_values.JxW(q);
- if (is_compressible)
+ if (scratch.material_model_outputs.is_compressible)
for (unsigned int i=0; i<dofs_per_cell; ++i)
data.local_pressure_shape_function_integrals(i) += scratch.phi_p[i] * scratch.finite_element_values.JxW(q);
}
@@ -1119,7 +1297,8 @@
?
update_gradients
:
- UpdateFlags(0)))),
+ UpdateFlags(0))),
+ parameters.n_compositional_fields),
internal::Assembly::CopyData::
StokesSystem<dim> (finite_element));
@@ -1164,7 +1343,14 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ compositional_fields.push_back(temp);
+ }
+
const unsigned int dofs_per_cell = scratch.finite_element_values.get_fe().dofs_per_cell;
const unsigned int n_q_points = scratch.finite_element_values.n_quadrature_points;
@@ -1203,16 +1389,17 @@
scratch.finite_element_values[pressure].get_function_values (old_old_solution,
scratch.old_old_pressure);
- scratch.finite_element_values[temperature].get_function_values(current_linearization_point,
- scratch.current_temperature_values);
scratch.finite_element_values[velocities].get_function_values(current_linearization_point,
scratch.current_velocity_values);
- scratch.finite_element_values[velocities].get_function_symmetric_gradients(current_linearization_point,
- scratch.current_strain_rates);
- scratch.finite_element_values[pressure].get_function_values(current_linearization_point,
- scratch.current_pressure_values);
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c) {
+ scratch.finite_element_values[compositional_fields[c]].get_function_values(old_solution,
+ scratch.old_composition_values[c]);
+ scratch.finite_element_values[compositional_fields[c]].get_function_values(old_old_solution,
+ scratch.old_old_composition_values[c]);
+ }
+
// TODO: Compute artificial viscosity once per timestep instead of each time
// temperature system is assembled (as this might happen more than once per
// timestep for iterative solvers)
@@ -1229,13 +1416,23 @@
scratch.old_old_strain_rates,
scratch.old_pressure,
scratch.old_old_pressure,
+ scratch.old_composition_values,
+ scratch.old_old_composition_values,
global_max_velocity,
global_T_range.second - global_T_range.first,
0.5 * (global_T_range.second + global_T_range.first),
global_entropy_variation,
scratch.finite_element_values.get_quadrature_points(),
- cell->diameter());
+ cell->diameter(),
+ 0); //index for temperature field
+ compute_material_model_input_values (current_linearization_point,
+ scratch.finite_element_values,
+ true,
+ scratch.material_model_inputs);
+
+ material_model->compute_parameters(scratch.material_model_inputs,scratch.material_model_outputs);
+
for (unsigned int q=0; q<n_q_points; ++q)
{
for (unsigned int k=0; k<dofs_per_cell; ++k)
@@ -1255,38 +1452,23 @@
:
scratch.old_temperature_values[q]);
- const double current_T = scratch.current_temperature_values[q];
- const SymmetricTensor<2,dim> current_strain_rate = scratch.current_strain_rates[q];
+ const double current_T = scratch.material_model_inputs.temperature[q];
+ const SymmetricTensor<2,dim> current_strain_rate = scratch.material_model_inputs.strain_rate[q];
const Tensor<1,dim> current_u = scratch.current_velocity_values[q];
- const double current_p = scratch.current_pressure_values[q];
+ const double current_p = scratch.material_model_inputs.pressure[q];
- const double alpha = material_model->thermal_expansion_coefficient(current_T,
- current_p,
- scratch.finite_element_values.quadrature_point(q));
- const double density = material_model->density(current_T,
- current_p,
- scratch.finite_element_values.quadrature_point(q));
- const double thermal_conductivity = material_model->thermal_conductivity(current_T,
- current_p,
- scratch.finite_element_values.quadrature_point(q));
- const double c_P = material_model->specific_heat (current_T,
- current_p,
- scratch.finite_element_values.quadrature_point(q));
+ const double alpha = scratch.material_model_outputs.thermal_expansion_coefficients[q];
+ const double density = scratch.material_model_outputs.densities[q];
+ const double thermal_conductivity = scratch.material_model_outputs.thermal_conductivities[q];
+ const double c_P = scratch.material_model_outputs.specific_heat[q];
+ const double viscosity = scratch.material_model_outputs.viscosities[q];
+ const bool is_compressible = scratch.material_model_outputs.is_compressible;
+ const double compressibility = (is_compressible
+ ?
+ scratch.material_model_outputs.compressibilities[q]
+ :
+ std::numeric_limits<double>::quiet_NaN() );
- const double viscosity = material_model->viscosity(current_T,
- current_p,
- current_strain_rate,
- scratch.finite_element_values.quadrature_point(q));
- const bool is_compressible = material_model->is_compressible ();
- const double compressibility
- = (is_compressible
- ?
- material_model->compressibility(scratch.old_temperature_values[q],
- scratch.old_pressure[q],
- scratch.finite_element_values.quadrature_point(q))
- :
- std::numeric_limits<double>::quiet_NaN() );
-
const Tensor<1,dim>
gravity = gravity_model->gravity_vector (scratch.finite_element_values.quadrature_point(q));
@@ -1412,7 +1594,8 @@
this,
std_cxx1x::_1),
internal::Assembly::Scratch::
- TemperatureSystem<dim> (finite_element, mapping, QGauss<dim>(parameters.temperature_degree+2)),
+ TemperatureSystem<dim> (finite_element, mapping, QGauss<dim>(parameters.temperature_degree+2),
+ parameters.n_compositional_fields),
internal::Assembly::CopyData::
TemperatureSystem<dim> (finite_element));
@@ -1486,24 +1669,9 @@
scratch.old_velocity_values);
scratch.finite_element_values[velocities].get_function_values (old_old_solution,
scratch.old_old_velocity_values);
- scratch.finite_element_values[velocities].get_function_symmetric_gradients (old_solution,
- scratch.old_strain_rates);
- scratch.finite_element_values[velocities].get_function_symmetric_gradients (old_old_solution,
- scratch.old_old_strain_rates);
- scratch.finite_element_values[pressure].get_function_values (old_solution,
- scratch.old_pressure);
- scratch.finite_element_values[pressure].get_function_values (old_old_solution,
- scratch.old_old_pressure);
-
- scratch.finite_element_values[composition].get_function_values(current_linearization_point,
- scratch.current_composition_values);
scratch.finite_element_values[velocities].get_function_values(current_linearization_point,
scratch.current_velocity_values);
- scratch.finite_element_values[velocities].get_function_symmetric_gradients(current_linearization_point,
- scratch.current_strain_rates);
- scratch.finite_element_values[pressure].get_function_values(current_linearization_point,
- scratch.current_pressure_values);
// TODO: Compute artificial viscosity once per timestep instead of each time
// temperature system is assembled (as this might happen more than once per
@@ -1517,16 +1685,19 @@
scratch.old_old_composition_laplacians,
scratch.old_velocity_values,
scratch.old_old_velocity_values,
- scratch.old_strain_rates,
- scratch.old_old_strain_rates,
- scratch.old_pressure,
- scratch.old_old_pressure,
+ std::vector<SymmetricTensor<2,dim> >(), //we do not need the strain rate for calculating the artificial viscosity
+ std::vector<SymmetricTensor<2,dim> >(), //strain rate
+ std::vector<double>(), //we also do not need the pressure
+ std::vector<double>(), //pressure
+ std::vector<std::vector<double> > (), //we also do not need the other compositional fields
+ std::vector<std::vector<double> > (), //other compositional fields
global_max_velocity,
global_C_range.second - global_C_range.first,
0.5 * (global_C_range.second + global_C_range.first),
global_entropy_variation,
- scratch.finite_element_values.get_quadrature_points(),
- cell->diameter());
+ std::vector<Point<dim, double> >(), //we also do not need the position for calculating the artificial viscosity
+ cell->diameter(),
+ composition_index+1); //index for the compositional field (0 is temperature)
for (unsigned int q=0; q<n_q_points; ++q)
{
@@ -1547,10 +1718,11 @@
:
scratch.old_composition_values[q]);
- const double current_C = scratch.current_composition_values[q];
const Tensor<1,dim> current_u = scratch.current_velocity_values[q];
- const double kappa = 1.e-12;
+ // the chemical diffusivity kappa is always smaller than the numerical diffusivity
+ // therefore we set it to zero
+ const double kappa = 0.0;
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
@@ -1634,7 +1806,8 @@
this,
std_cxx1x::_1),
internal::Assembly::Scratch::
- CompositionSystem<dim> (finite_element, mapping, QGauss<dim>(parameters.composition_degree+2)),
+ CompositionSystem<dim> (finite_element, mapping, QGauss<dim>(parameters.composition_degree+2),
+ parameters.n_compositional_fields),
internal::Assembly::CopyData::
CompositionSystem<dim> (finite_element));
Modified: trunk/aspect/source/simulator/core.cc
===================================================================
--- trunk/aspect/source/simulator/core.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/simulator/core.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -103,9 +103,7 @@
*boundary_temperature,
*adiabatic_conditions)),
compositional_initial_conditions (CompositionalInitialConditions::create_initial_conditions (prm,
- *geometry_model,
- *boundary_temperature,
- *adiabatic_conditions)),
+ *geometry_model)),
time (std::numeric_limits<double>::quiet_NaN()),
time_step (0),
@@ -218,8 +216,10 @@
adiabatic_conditions.reset (new AdiabaticConditions<dim>(*geometry_model,
*gravity_model,
*material_model,
+ *compositional_initial_conditions,
parameters.surface_pressure,
- parameters.adiabatic_surface_temperature));
+ parameters.adiabatic_surface_temperature,
+ parameters.n_compositional_fields));
for (std::map<types::boundary_id_t,std::string>::const_iterator
p = parameters.prescribed_velocity_boundary_indicators.begin();
p != parameters.prescribed_velocity_boundary_indicators.end();
@@ -571,10 +571,10 @@
n_T = system_dofs_per_block[2];
unsigned int n_C_sum = 0;
std::vector<unsigned int> n_C (parameters.n_compositional_fields+1);
- for (unsigned int i=0; i<parameters.n_compositional_fields; ++i)
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
{
- n_C[i] = system_dofs_per_block[i+3];
- n_C_sum += n_C[i];
+ n_C[c] = system_dofs_per_block[c+3];
+ n_C_sum += n_C[c];
}
@@ -602,8 +602,8 @@
<< n_u + n_p + n_T + n_C_sum
<< " (" << n_u << '+' << n_p << '+'<< n_T;
- for (unsigned int i=0; i<parameters.n_compositional_fields; ++i)
- pcout << '+' << n_C[i];
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+ pcout << '+' << n_C[c];
pcout <<')'
<< std::endl
@@ -631,13 +631,13 @@
{
unsigned int n_C_so_far = 0;
- for (unsigned int i=0; i<parameters.n_compositional_fields; ++i)
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
{
system_partitioning.push_back(system_index_set.get_view(n_u+n_p+n_T+n_C_so_far,
- n_u+n_p+n_T+n_C_so_far+n_C[i]));
+ n_u+n_p+n_T+n_C_so_far+n_C[c]));
system_relevant_partitioning.push_back(system_relevant_set.get_view(n_u+n_p+n_T+n_C_so_far,
- n_u+n_p+n_T+n_C_so_far+n_C[i]));
- n_C_so_far += n_C[i];
+ n_u+n_p+n_T+n_C_so_far+n_C[c]));
+ n_C_so_far += n_C[c];
}
}
}
@@ -787,14 +787,24 @@
Vector<float> estimated_error_per_cell_rho (triangulation.n_active_cells());
Vector<float> estimated_error_per_cell_T (triangulation.n_active_cells());
Vector<float> estimated_error_per_cell_u (triangulation.n_active_cells());
+ std::vector<Vector<float> > estimated_error_per_cell_C (parameters.n_compositional_fields,
+ Vector<float> (triangulation.n_active_cells()));
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> composition;
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ composition.push_back(temp);
+ }
+
//Velocity|Temperature|Normalized density and temperature|Weighted density and temperature|Density c_p temperature
// compute density error
- if (parameters.refinement_strategy != "Temperature" && parameters.refinement_strategy != "Velocity")
+ if (parameters.refinement_strategy != "Temperature" && parameters.refinement_strategy != "Velocity"
+ && parameters.refinement_strategy != "Composition")
{
bool lookup_rho_c_p_T = (parameters.refinement_strategy == "Density c_p temperature");
@@ -818,7 +828,14 @@
std::vector<double> pressure_values(quadrature.size());
std::vector<double> temperature_values(quadrature.size());
+ // the values of the compositional fields are stored as blockvectors for each field
+ // we have to extract them in this structure
+ std::vector<std::vector<double> > prelim_composition_values (parameters.n_compositional_fields,
+ std::vector<double> (quadrature.size()));
+ std::vector<std::vector<double> > composition_values (quadrature.size(),
+ std::vector<double> (parameters.n_compositional_fields));
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
@@ -830,6 +847,14 @@
pressure_values);
fe_values[temperature].get_function_values (solution,
temperature_values);
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ fe_values[composition[c]].get_function_values (solution,
+ prelim_composition_values[c]);
+ // then we copy these values to exchange the inner and outer vector, because for the material
+ // model we need a vector with values of all the compositional fields for every quadrature point
+ for(unsigned int q=0;q<quadrature.size();++q)
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ composition_values[q][c] = prelim_composition_values[c][q];
cell->get_dof_indices (local_dof_indices);
@@ -845,12 +870,14 @@
vec_distributed(local_dof_indices[system_local_dof])
= material_model->density( temperature_values[i],
pressure_values[i],
+ composition_values[i],
fe_values.quadrature_point(i))
* ((lookup_rho_c_p_T)
?
(temperature_values[i]
* material_model->specific_heat(temperature_values[i],
pressure_values[i],
+ composition_values[i],
fe_values.quadrature_point(i)))
:
1.0);
@@ -914,6 +941,30 @@
estimated_error_per_cell_T = 0;
}
+ // compute the errors for composition solution
+ if (parameters.refinement_strategy == "Composition")
+ {
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ std::vector<bool> composition_component (dim+2+parameters.n_compositional_fields, false);
+ composition_component[dim+2+c] = true;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ QGauss<dim-1>(parameters.composition_degree+1),
+ typename FunctionMap<dim>::type(),
+ solution,
+ estimated_error_per_cell_C[c],
+ composition_component,
+ 0,
+ 0,
+ triangulation.locally_owned_subdomain());
+ }
+ }
+ else
+ {
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ estimated_error_per_cell_C[c] = 0;
+ }
+
// compute the errors for the stokes solution
if (parameters.refinement_strategy == "Velocity")
{
@@ -983,6 +1034,14 @@
for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
estimated_error_per_cell(i) = estimated_error_per_cell_rho(i);
}
+ else if (parameters.refinement_strategy == "Composition")
+ {
+ for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+ estimated_error_per_cell(i) = 0;
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+ for (unsigned int i=0; i<estimated_error_per_cell.size(); ++i)
+ estimated_error_per_cell(i) += estimated_error_per_cell_C[c](i);
+ }
else
AssertThrow(false, ExcNotImplemented());
}
@@ -1088,12 +1147,12 @@
current_linearization_point.block(2) = solution.block(2);
- for (unsigned int n=0; n<parameters.n_compositional_fields; ++n)
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
{
- assemble_composition_system (n);
- build_composition_preconditioner(n);
- solve_temperature_or_composition(1+n); // this is correct, 0 would be temperature
- current_linearization_point.block(3+n) = solution.block(3+n);
+ assemble_composition_system (c);
+ build_composition_preconditioner(c);
+ solve_temperature_or_composition(1+c); // this is correct, 0 would be temperature
+ current_linearization_point.block(3+c) = solution.block(3+c);
}
assemble_stokes_system();
@@ -1124,13 +1183,13 @@
rebuild_stokes_matrix = true;
std::vector<double> composition_residual (parameters.n_compositional_fields,0);
- for (unsigned int n=0; n<parameters.n_compositional_fields; ++n)
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
{
- assemble_composition_system (n);
- build_composition_preconditioner(n);
- composition_residual[n]
- = solve_temperature_or_composition(1+n); // 1+n is correct, because 0 is for temperature
- current_linearization_point.block(3+n) = solution.block(3+n);
+ assemble_composition_system (c);
+ build_composition_preconditioner(c);
+ composition_residual[c]
+ = solve_temperature_or_composition(1+c); // 1+n is correct, because 0 is for temperature
+ current_linearization_point.block(3+c) = solution.block(3+c);
}
assemble_stokes_system();
@@ -1144,8 +1203,8 @@
pcout << " Nonlinear residuals: " << temperature_residual
<< ", " << stokes_residual;
- for (unsigned int n=0; n<parameters.n_compositional_fields; ++n)
- pcout << ", " << composition_residual[n];
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+ pcout << ", " << composition_residual[c];
pcout << std::endl
<< std::endl;
@@ -1153,16 +1212,16 @@
if (iteration == 0)
{
initial_temperature_residual = temperature_residual;
- for (unsigned int n=0; n<parameters.n_compositional_fields; ++n)
- initial_composition_residual[n] = composition_residual[n];
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+ initial_composition_residual[c] = composition_residual[c];
initial_stokes_residual = stokes_residual;
}
else
{
//TODO: make this a parameter in the input file
double max = 0.0;
- for (unsigned int n=0; n<parameters.n_compositional_fields; ++n)
- max = std::max(composition_residual[n]/initial_composition_residual[n],max);
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+ max = std::max(composition_residual[c]/initial_composition_residual[c],max);
if (std::max(std::max(stokes_residual/initial_stokes_residual,
temperature_residual/initial_temperature_residual),
max) < 1e-3)
@@ -1183,10 +1242,10 @@
assemble_temperature_system ();
solve_temperature_or_composition(0);
- for (unsigned int n=0; n<parameters.n_compositional_fields; ++n)
+ for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
{
- assemble_composition_system (n);
- solve_temperature_or_composition(1+n); // 1+n is correct, because 0 is the temperature
+ assemble_composition_system (c);
+ solve_temperature_or_composition(1+c); // 1+n is correct, because 0 is the temperature
}
// ...and then iterate the solution
Modified: trunk/aspect/source/simulator/helper_functions.cc
===================================================================
--- trunk/aspect/source/simulator/helper_functions.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/simulator/helper_functions.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -254,7 +254,16 @@
}
+ template <int dim>
+ void Simulator<dim>::extract_composition_values_at_q_point (const std::vector<std::vector<double> > &composition_values,
+ const unsigned int q,
+ std::vector<double> &composition_values_at_q_point) const
+ {
+ for(unsigned int k=0; k < composition_values_at_q_point.size(); ++k)
+ composition_values_at_q_point[k] = composition_values[k][q];
+ }
+
template <int dim>
std::pair<double,double>
Simulator<dim>::get_extrapolated_temperature_or_composition_range (const unsigned int index) const
@@ -616,10 +625,19 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ compositional_fields.push_back(temp);
+ }
+
std::vector<SymmetricTensor<2,dim> > strain_rates(n_q_points);
std::vector<double> pressure_values(n_q_points);
std::vector<double> temperature_values(n_q_points);
+ std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,std::vector<double> (n_q_points));
+ std::vector<double> composition_values_at_q_point (parameters.n_compositional_fields);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
@@ -637,6 +655,9 @@
temperature_values);
fe_values[velocities].get_function_symmetric_gradients (this->solution,
strain_rates);
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ fe_values[compositional_fields[c]].get_function_values(this->solution,
+ composition_values[c]);
for (unsigned int q = 0; q < n_q_points; ++q)
{
@@ -644,8 +665,15 @@
const unsigned int idx = static_cast<unsigned int>((depth*num_slices)/max_depth);
Assert(idx<num_slices, ExcInternalError());
+ extract_composition_values_at_q_point (composition_values,
+ q,
+ composition_values_at_q_point);
+
+ // TODO: we should use compute_parameters() instead. This should be done
+ // together with merging all the *_average() functions.
const double viscosity = material_model->viscosity(temperature_values[q],
pressure_values[q],
+ composition_values_at_q_point,
strain_rates[q],
fe_values.quadrature_point(q));
++counts[idx];
@@ -799,8 +827,17 @@
update_values | update_quadrature_points);
const FEValuesExtractors::Scalar pressure (dim);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ compositional_fields.push_back(temp);
+ }
+
std::vector<double> pressure_values(n_q_points);
+ std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,std::vector<double> (n_q_points));
+ std::vector<double> composition_values_at_q_point (parameters.n_compositional_fields);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
@@ -814,14 +851,22 @@
fe_values.reinit (cell);
fe_values[pressure].get_function_values (this->solution,
pressure_values);
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ fe_values[compositional_fields[c]].get_function_values(this->solution,
+ composition_values[c]);
for (unsigned int q = 0; q < n_q_points; ++q)
{
const double depth = geometry_model->depth(fe_values.quadrature_point(q));
const unsigned int idx = static_cast<unsigned int>((depth*num_slices)/max_depth);
Assert(idx<num_slices, ExcInternalError());
+ extract_composition_values_at_q_point (composition_values,
+ q,
+ composition_values_at_q_point);
+
const double Vs_depth_average = material_model->seismic_Vs(average_temperature[idx],
pressure_values[q],
+ composition_values_at_q_point,
fe_values.quadrature_point(q));
++counts[idx];
values[idx] += Vs_depth_average;
@@ -863,8 +908,17 @@
update_quadrature_points );
const FEValuesExtractors::Scalar pressure (dim);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ compositional_fields.push_back(temp);
+ }
+
std::vector<double> pressure_values(n_q_points);
+ std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,std::vector<double> (n_q_points));
+ std::vector<double> composition_values_at_q_point (parameters.n_compositional_fields);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
@@ -878,14 +932,22 @@
fe_values.reinit (cell);
fe_values[pressure].get_function_values (this->solution,
pressure_values);
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ fe_values[compositional_fields[c]].get_function_values(this->solution,
+ composition_values[c]);
for (unsigned int q = 0; q < n_q_points; ++q)
{
const double depth = geometry_model->depth(fe_values.quadrature_point(q));
const unsigned int idx = static_cast<unsigned int>((depth*num_slices)/max_depth);
Assert(idx<num_slices, ExcInternalError());
+ extract_composition_values_at_q_point (composition_values,
+ q,
+ composition_values_at_q_point);
+
const double Vp_depth_average = material_model->seismic_Vp(average_temperature[idx],
pressure_values[q],
+ composition_values_at_q_point,
fe_values.quadrature_point(q));
++counts[idx];
values[idx] += Vp_depth_average;
@@ -927,9 +989,18 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ compositional_fields.push_back(temp);
+ }
+
std::vector<double> pressure_values(n_q_points);
std::vector<double> temperature_values(n_q_points);
+ std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,std::vector<double> (n_q_points));
+ std::vector<double> composition_values_at_q_point (parameters.n_compositional_fields);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
@@ -944,9 +1015,18 @@
pressure_values);
fe_values[temperature].get_function_values (this->solution,
temperature_values);
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ fe_values[compositional_fields[c]].get_function_values(this->solution,
+ composition_values[c]);
+
+ extract_composition_values_at_q_point (composition_values,
+ 0,
+ composition_values_at_q_point);
+
const double Vs = material_model->seismic_Vs(temperature_values[0],
pressure_values[0],
+ composition_values_at_q_point,
fe_values.quadrature_point(0));
const double depth = geometry_model->depth(fe_values.quadrature_point(0));
const unsigned int idx = static_cast<unsigned int>((depth*num_slices)/max_depth);
@@ -984,9 +1064,18 @@
const FEValuesExtractors::Vector velocities (0);
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+c);
+ compositional_fields.push_back(temp);
+ }
+
std::vector<double> pressure_values(n_q_points);
std::vector<double> temperature_values(n_q_points);
+ std::vector<std::vector<double> > composition_values (parameters.n_compositional_fields,std::vector<double> (n_q_points));
+ std::vector<double> composition_values_at_q_point (parameters.n_compositional_fields);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
@@ -1002,9 +1091,17 @@
pressure_values);
fe_values[temperature].get_function_values (this->solution,
temperature_values);
+ for(unsigned int c=0;c<parameters.n_compositional_fields;++c)
+ fe_values[compositional_fields[c]].get_function_values(this->solution,
+ composition_values[c]);
+ extract_composition_values_at_q_point (composition_values,
+ 0,
+ composition_values_at_q_point);
+
const double Vp = material_model->seismic_Vp(temperature_values[0],
pressure_values[0],
+ composition_values_at_q_point,
fe_values.quadrature_point(0));
const double depth = geometry_model->depth(fe_values.quadrature_point(0));
const unsigned int idx = static_cast<unsigned int>((depth*num_slices)/max_depth);
@@ -1021,6 +1118,9 @@
template void Simulator<dim>::normalize_pressure(LinearAlgebra::BlockVector &vector); \
template double Simulator<dim>::get_maximal_velocity (const LinearAlgebra::BlockVector &solution) const; \
template std::pair<double,double> Simulator<dim>::get_extrapolated_temperature_or_composition_range (const unsigned int index) const; \
+ template void Simulator<dim>::extract_composition_values_at_q_point (const std::vector<std::vector<double> > &composition_values, \
+ const unsigned int q, \
+ std::vector<double> &composition_values_at_q_point) const; \
template double Simulator<dim>::compute_time_step () const; \
template void Simulator<dim>::make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector); \
template void Simulator<dim>::compute_depth_average_field(const unsigned int index, std::vector<double> &values) const; \
Modified: trunk/aspect/source/simulator/parameters.cc
===================================================================
--- trunk/aspect/source/simulator/parameters.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/simulator/parameters.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -299,7 +299,8 @@
"Velocity|"
"Normalized density and temperature|"
"Weighted density and temperature|"
- "Density c_p temperature"),
+ "Density c_p temperature|"
+ "Composition"),
"The method used to determine which cells to refine and which "
"to coarsen.");
prm.declare_entry ("Additional refinement times", "",
Modified: trunk/aspect/source/simulator/simulator_access.cc
===================================================================
--- trunk/aspect/source/simulator/simulator_access.cc 2012-10-19 15:09:58 UTC (rev 1296)
+++ trunk/aspect/source/simulator/simulator_access.cc 2012-10-19 15:30:58 UTC (rev 1297)
@@ -287,6 +287,18 @@
{
return *simulator->adiabatic_conditions.get();
}
+
+ template <int dim>
+ void
+ SimulatorAccess<dim>::get_composition_values_at_q_point (const std::vector<std::vector<double>> &composition_values,
+ const unsigned int q,
+ std::vector<double> &composition_values_at_q_point) const
+ {
+ simulator->extract_composition_values_at_q_point (composition_values,
+ q,
+ composition_values_at_q_point);
+ }
+
}
More information about the CIG-COMMITS
mailing list