[cig-commits] r1279 - in branches/active_compositions: 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
Tue Oct 16 11:05:49 PDT 2012
Author: dannberg
Date: 2012-10-16 12:05:48 -0600 (Tue, 16 Oct 2012)
New Revision: 1279
Modified:
branches/active_compositions/include/aspect/adiabatic_conditions.h
branches/active_compositions/include/aspect/compositional_initial_conditions/interface.h
branches/active_compositions/include/aspect/material_model/duretz_et_al.h
branches/active_compositions/include/aspect/material_model/interface.h
branches/active_compositions/include/aspect/material_model/simple.h
branches/active_compositions/include/aspect/material_model/steinberger.h
branches/active_compositions/include/aspect/material_model/table.h
branches/active_compositions/include/aspect/material_model/tan_gurnis.h
branches/active_compositions/include/aspect/simulator.h
branches/active_compositions/source/adiabatic_conditions.cc
branches/active_compositions/source/compositional_initial_conditions/interface.cc
branches/active_compositions/source/material_model/duretz_et_al.cc
branches/active_compositions/source/material_model/interface.cc
branches/active_compositions/source/material_model/simple.cc
branches/active_compositions/source/material_model/steinberger.cc
branches/active_compositions/source/material_model/table.cc
branches/active_compositions/source/material_model/tan_gurnis.cc
branches/active_compositions/source/postprocess/heat_flux_statistics.cc
branches/active_compositions/source/postprocess/table_heat_flux_statistics.cc
branches/active_compositions/source/postprocess/table_velocity_statistics.cc
branches/active_compositions/source/postprocess/velocity_statistics.cc
branches/active_compositions/source/postprocess/visualization/density.cc
branches/active_compositions/source/postprocess/visualization/seismic_vp.cc
branches/active_compositions/source/postprocess/visualization/seismic_vs.cc
branches/active_compositions/source/postprocess/visualization/specific_heat.cc
branches/active_compositions/source/postprocess/visualization/thermal_expansivity.cc
branches/active_compositions/source/postprocess/visualization/thermodynamic_phase.cc
branches/active_compositions/source/simulator/core.cc
branches/active_compositions/source/simulator/helper_functions.cc
Log:
Make all properties in all of the material models dependent on the compositional fields = active composition.
Modified: branches/active_compositions/include/aspect/adiabatic_conditions.h
===================================================================
--- branches/active_compositions/include/aspect/adiabatic_conditions.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/adiabatic_conditions.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -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: branches/active_compositions/include/aspect/compositional_initial_conditions/interface.h
===================================================================
--- branches/active_compositions/include/aspect/compositional_initial_conditions/interface.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/compositional_initial_conditions/interface.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -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>
@@ -65,8 +64,7 @@
*/
void
initialize (const GeometryModel::Interface<dim> &geometry_model,
- const BoundaryTemperature::Interface<dim> &boundary_temperature,
- const AdiabaticConditions<dim> &adiabatic_conditions);
+ const BoundaryTemperature::Interface<dim> &boundary_temperature);
/**
* Return the initial temperature as a function of position.
@@ -107,10 +105,6 @@
*/
const BoundaryTemperature::Interface<dim> *boundary_temperature;
- /**
- * Pointer to an object that describes adiabatic conditions.
- */
- const AdiabaticConditions<dim> *adiabatic_conditions;
};
@@ -151,8 +145,7 @@
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 BoundaryTemperature::Interface<dim> &boundary_temperature);
/**
Modified: branches/active_compositions/include/aspect/material_model/duretz_et_al.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/duretz_et_al.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/material_model/duretz_et_al.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -82,22 +82,27 @@
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;
/**
* @}
@@ -251,22 +256,27 @@
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;
/**
* @}
@@ -373,22 +383,27 @@
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: branches/active_compositions/include/aspect/material_model/interface.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/interface.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/material_model/interface.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -118,6 +118,7 @@
*/
virtual double density (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const = 0;
/**
@@ -127,6 +128,7 @@
*/
virtual double compressibility (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const = 0;
/**
@@ -135,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;
/**
@@ -151,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;
/**
@@ -169,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;
/**
* @}
@@ -247,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;
@@ -261,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;
@@ -275,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;
@@ -289,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;
@@ -303,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;
/**
@@ -362,6 +372,7 @@
double
seismic_Vp (const double temperature,
const double pressure,
+ const std::vector<double> &compositional_fields,
const Point<dim> &position) const;
/**
@@ -379,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
@@ -393,7 +405,8 @@
virtual
unsigned int
thermodynamic_phase (const double temperature,
- const double pressure) const;
+ const double pressure,
+ const std::vector<double> &compositional_fields) const;
/**
* @}
*/
Modified: branches/active_compositions/include/aspect/material_model/simple.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/simple.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/material_model/simple.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -59,26 +59,28 @@
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;
-
- virtual void compute_parameters(typename Interface<dim>::MaterialModelInputs &in,
- typename Interface<dim>::MaterialModelOutputs &out);
/**
* @}
*/
Modified: branches/active_compositions/include/aspect/material_model/steinberger.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/steinberger.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/material_model/steinberger.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -62,30 +62,37 @@
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: branches/active_compositions/include/aspect/material_model/table.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/table.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/material_model/table.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -60,18 +60,22 @@
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;
/**
* @}
@@ -182,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;
/**
@@ -209,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
@@ -219,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
@@ -231,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: branches/active_compositions/include/aspect/material_model/tan_gurnis.h
===================================================================
--- branches/active_compositions/include/aspect/material_model/tan_gurnis.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/material_model/tan_gurnis.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -57,22 +57,27 @@
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: branches/active_compositions/include/aspect/simulator.h
===================================================================
--- branches/active_compositions/include/aspect/simulator.h 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/include/aspect/simulator.h 2012-10-16 18:05:48 UTC (rev 1279)
@@ -1364,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: branches/active_compositions/source/adiabatic_conditions.cc
===================================================================
--- branches/active_compositions/source/adiabatic_conditions.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/adiabatic_conditions.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -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,20 @@
const Point<dim> representative_point = geometry_model.representative_point (z);
+ std::vector<double> initial_composition(n_compositional_fields);
+ for (unsigned int k=0;k<n_compositional_fields;++k)
+ initial_composition[k] = compositional_initial_conditions.initial_composition(representative_point, k);
+
// 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]
Modified: branches/active_compositions/source/compositional_initial_conditions/interface.cc
===================================================================
--- branches/active_compositions/source/compositional_initial_conditions/interface.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/compositional_initial_conditions/interface.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -41,12 +41,11 @@
template <int dim>
void
Interface<dim>::initialize (const GeometryModel::Interface<dim> &geometry_model_,
- const BoundaryTemperature::Interface<dim> &boundary_temperature_,
- const AdiabaticConditions<dim> &adiabatic_conditions_)
+ const BoundaryTemperature::Interface<dim> &boundary_temperature_)
+ // TODO check why we need boundary temperature. is this boundary composition?
{
geometry_model = &geometry_model_;
boundary_temperature = &boundary_temperature_;
- adiabatic_conditions = &adiabatic_conditions_;
}
@@ -95,8 +94,7 @@
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 BoundaryTemperature::Interface<dim> &boundary_temperature)
{
// we need to get at the number of compositional fields here to
// know whether we need to generate something at all.
@@ -117,8 +115,7 @@
Interface<dim> *plugin = std_cxx1x::get<dim>(registered_plugins).create_plugin (model_name, prm);
plugin->initialize (geometry_model,
- boundary_temperature,
- adiabatic_conditions);
+ boundary_temperature);
return plugin;
}
}
@@ -183,8 +180,7 @@
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 BoundaryTemperature::Interface<dim> &boundary_temperature);
ASPECT_INSTANTIATE(INSTANTIATE)
}
Modified: branches/active_compositions/source/material_model/duretz_et_al.cc
===================================================================
--- branches/active_compositions/source/material_model/duretz_et_al.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/material_model/duretz_et_al.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -76,6 +76,7 @@
SolCx<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -94,6 +95,7 @@
SolCx<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -112,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
@@ -125,6 +128,7 @@
SolCx<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -136,6 +140,7 @@
SolCx<dim>::
compressibility (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0.0;
@@ -297,6 +302,7 @@
SolKz<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -315,6 +321,7 @@
SolKz<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -333,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
@@ -345,6 +353,7 @@
SolKz<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -356,6 +365,7 @@
SolKz<dim>::
compressibility (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0.0;
@@ -460,6 +470,7 @@
Inclusion<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -478,6 +489,7 @@
Inclusion<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -496,6 +508,7 @@
Inclusion<dim>::
density (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &p) const
{
return 0;
@@ -507,6 +520,7 @@
Inclusion<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0;
@@ -518,6 +532,7 @@
Inclusion<dim>::
compressibility (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0.0;
Modified: branches/active_compositions/source/material_model/interface.cc
===================================================================
--- branches/active_compositions/source/material_model/interface.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/material_model/interface.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -48,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
{
@@ -62,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
{
@@ -75,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
{
@@ -88,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
{
@@ -101,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
{
@@ -185,6 +190,7 @@
Interface<dim>::
seismic_Vp (double dummy1,
double dummy2,
+ const std::vector<double> &, /*composition*/
const Point<dim> &dummy3) const
{
return -1.0;
@@ -196,6 +202,7 @@
Interface<dim>::
seismic_Vs (double dummy1,
double dummy2,
+ const std::vector<double> &, /*composition*/
const Point<dim> &dummy3) const
{
return -1.0;
@@ -206,7 +213,8 @@
unsigned int
Interface<dim>::
thermodynamic_phase (double dummy1,
- double dummy2) const
+ double dummy2,
+ const std::vector<double> & /*composition*/) const
{
return 0;
}
@@ -216,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>
@@ -276,14 +285,14 @@
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.position[i]);
- out.thermal_expansion_coefficients[i] = thermal_expansion_coefficient (in.temperature[i], in.pressure[i], in.position[i]);
- out.seismic_Vp[i] = seismic_Vp (in.temperature[i], in.pressure[i], in.position[i]);
- out.seismic_Vs[i] = seismic_Vs (in.temperature[i], in.pressure[i], in.position[i]);
- out.specific_heat[i] = specific_heat (in.temperature[i], in.pressure[i], in.position[i]);
- out.thermal_conductivities[i] = thermal_conductivity (in.temperature[i], in.pressure[i], in.position[i]);
- out.compressibilities[i] = compressibility (in.temperature[i], in.pressure[i], in.position[i]);
- out.thermodynamic_phases[i] = thermodynamic_phase (in.temperature[i], in.pressure[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.seismic_Vp[i] = seismic_Vp (in.temperature[i], in.pressure[i], in.composition[i], in.position[i]);
+ out.seismic_Vs[i] = seismic_Vs (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.thermodynamic_phases[i] = thermodynamic_phase (in.temperature[i], in.pressure[i], in.composition[i]);
out.is_compressible = is_compressible();
}
}
Modified: branches/active_compositions/source/material_model/simple.cc
===================================================================
--- branches/active_compositions/source/material_model/simple.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/material_model/simple.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -40,7 +40,11 @@
const Point<dim> &) const
{
// return eta;
- return (4.0*composition[0]+1) * eta;
+ return (this->n_compositional_fields()>0
+ ?
+ (4.0*composition[0]+1) * eta
+ :
+ eta);
}
@@ -73,6 +77,7 @@
Simple<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return reference_specific_heat;
@@ -91,6 +96,7 @@
Simple<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return k_value;
@@ -108,11 +114,16 @@
double
Simple<dim>::
density (const double temperature,
- const double comp,
+ 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)));
}
@@ -121,6 +132,7 @@
Simple<dim>::
thermal_expansion_coefficient (const double temperature,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return thermal_alpha;
@@ -132,37 +144,13 @@
Simple<dim>::
compressibility (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 0.0;
}
template <int dim>
- void
- Simple<dim>::compute_parameters(struct Interface<dim>::MaterialModelInputs &in, struct Interface<dim>::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] = (this->n_compositional_fields()>0)
- ?
- (100*in.composition[i][0]) + reference_rho *(1 - thermal_alpha * (in.temperature[i] - reference_T))
- :
- density (in.temperature[i], in.pressure[i], in.position[i]);
- out.thermal_expansion_coefficients[i] = thermal_expansion_coefficient (in.temperature[i], in.pressure[i], in.position[i]);
- out.seismic_Vp[i] = seismic_Vp (in.temperature[i], in.pressure[i], in.position[i]);
- out.seismic_Vs[i] = seismic_Vs (in.temperature[i], in.pressure[i], in.position[i]);
- out.specific_heat[i] = specific_heat (in.temperature[i], in.pressure[i], in.position[i]);
- out.thermal_conductivities[i] = thermal_conductivity (in.temperature[i], in.pressure[i], in.position[i]);
- out.compressibilities[i] = compressibility (in.temperature[i], in.pressure[i], in.position[i]);
- out.thermodynamic_phases[i] = thermodynamic_phase (in.temperature[i], in.pressure[i]);
- out.is_compressible = is_compressible();
- }
- }
-
-
-
- template <int dim>
bool
Simple<dim>::
viscosity_depends_on (const NonlinearDependence::Dependence) const
Modified: branches/active_compositions/source/material_model/steinberger.cc
===================================================================
--- branches/active_compositions/source/material_model/steinberger.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/material_model/steinberger.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -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,7 +411,7 @@
Steinberger<dim>::
viscosity (const double temperature,
const double,
- const std::vector<double> &, /*composition*/
+ const std::vector<double> &compositional_fields,
const SymmetricTensor<2,dim> &,
const Point<dim> &position) const
{
@@ -464,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;
@@ -484,6 +505,7 @@
Steinberger<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 4.7;
@@ -495,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;
@@ -506,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;
@@ -515,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);
}
}
@@ -525,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;
@@ -536,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;
@@ -548,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: branches/active_compositions/source/material_model/table.cc
===================================================================
--- branches/active_compositions/source/material_model/table.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/material_model/table.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -468,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");
@@ -479,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 */ //??
@@ -502,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
@@ -521,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");
@@ -531,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);
@@ -542,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);
@@ -553,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;
@@ -568,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: branches/active_compositions/source/material_model/tan_gurnis.cc
===================================================================
--- branches/active_compositions/source/material_model/tan_gurnis.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/material_model/tan_gurnis.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -88,6 +88,7 @@
TanGurnis<dim>::
specific_heat (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 1250;
@@ -106,6 +107,7 @@
TanGurnis<dim>::
thermal_conductivity (const double,
const double,
+ const std::vector<double> &, /*composition*/
const Point<dim> &) const
{
return 2e-5;
@@ -124,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);
@@ -137,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;
@@ -148,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: branches/active_compositions/source/postprocess/heat_flux_statistics.cc
===================================================================
--- branches/active_compositions/source/postprocess/heat_flux_statistics.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/heat_flux_statistics.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -53,10 +53,19 @@
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int q=0;q<this->n_compositional_fields();++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ 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 i=0;i<this->n_compositional_fields();++i)
+ fe_face_values[compositional_fields[i]].get_function_values(this->get_solution(),
+ composition_values[i]);
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: branches/active_compositions/source/postprocess/table_heat_flux_statistics.cc
===================================================================
--- branches/active_compositions/source/postprocess/table_heat_flux_statistics.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/table_heat_flux_statistics.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -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 q=0;q<this->n_compositional_fields();++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ 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 i=0;i<this->n_compositional_fields();++i)
+ fe_face_values[compositional_fields[i]].get_function_values(this->get_solution(),
+ composition_values[i]);
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: branches/active_compositions/source/postprocess/table_velocity_statistics.cc
===================================================================
--- branches/active_compositions/source/postprocess/table_velocity_statistics.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/table_velocity_statistics.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -165,6 +165,7 @@
if (this->get_time() == 0e0)
{
double dT = this->get_boundary_temperature().maximal_temperature() - this->get_boundary_temperature().minimal_temperature();
+ 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 +202,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)
<< 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())
<< std::endl;
this->get_pcout()<< std::endl;
}
Modified: branches/active_compositions/source/postprocess/velocity_statistics.cc
===================================================================
--- branches/active_compositions/source/postprocess/velocity_statistics.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/velocity_statistics.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -148,6 +148,7 @@
const double h = this->get_geometry_model().maximal_depth();
const double dT = this->get_boundary_temperature().maximal_temperature() - this->get_boundary_temperature().minimal_temperature();
+ 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 +185,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)
<< 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())
<< std::endl;
this->get_pcout()<< std::endl;
}
Modified: branches/active_compositions/source/postprocess/visualization/density.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/density.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/visualization/density.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -68,6 +68,7 @@
computed_quantities[q](0) = this->get_material_model().density(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: branches/active_compositions/source/postprocess/visualization/seismic_vp.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/seismic_vp.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/visualization/seismic_vp.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -67,6 +67,7 @@
computed_quantities[q](0) = this->get_material_model().seismic_Vp(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: branches/active_compositions/source/postprocess/visualization/seismic_vs.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/seismic_vs.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/visualization/seismic_vs.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -67,6 +67,7 @@
computed_quantities[q](0) = this->get_material_model().seismic_Vs(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: branches/active_compositions/source/postprocess/visualization/specific_heat.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/specific_heat.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/visualization/specific_heat.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -67,6 +67,7 @@
computed_quantities[q](0) = this->get_material_model().specific_heat(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: branches/active_compositions/source/postprocess/visualization/thermal_expansivity.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/thermal_expansivity.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/visualization/thermal_expansivity.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -68,6 +68,7 @@
computed_quantities[q](0) = this->get_material_model().thermal_expansion_coefficient(temperature,
pressure,
+ composition,
evaluation_points[q]);
}
}
Modified: branches/active_compositions/source/postprocess/visualization/thermodynamic_phase.cc
===================================================================
--- branches/active_compositions/source/postprocess/visualization/thermodynamic_phase.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/postprocess/visualization/thermodynamic_phase.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -66,7 +66,8 @@
composition[i] = uh[q][dim+2+i];
computed_quantities[q](0) = this->get_material_model().thermodynamic_phase(temperature,
- pressure);
+ pressure,
+ composition);
}
}
}
Modified: branches/active_compositions/source/simulator/core.cc
===================================================================
--- branches/active_compositions/source/simulator/core.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/simulator/core.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -104,8 +104,7 @@
*adiabatic_conditions)),
compositional_initial_conditions (CompositionalInitialConditions::create_initial_conditions (prm,
*geometry_model,
- *boundary_temperature,
- *adiabatic_conditions)),
+ *boundary_temperature)),
time (std::numeric_limits<double>::quiet_NaN()),
time_step (0),
@@ -218,8 +217,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();
@@ -790,7 +791,14 @@
const FEValuesExtractors::Scalar pressure (dim);
const FEValuesExtractors::Scalar temperature (dim+1);
+ std::vector<FEValuesExtractors::Scalar> composition;
+ for (unsigned int q=0;q<parameters.n_compositional_fields;++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ composition.push_back(temp);
+ }
+
//Velocity|Temperature|Normalized density and temperature|Weighted density and temperature|Density c_p temperature
// compute density error
@@ -818,7 +826,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 +845,14 @@
pressure_values);
fe_values[temperature].get_function_values (solution,
temperature_values);
+ for(unsigned int i=0;i<parameters.n_compositional_fields;++i)
+ fe_values[composition[i]].get_function_values (solution,
+ prelim_composition_values[i]);
+ // 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 i=0;i<parameters.n_compositional_fields;++i)
+ composition_values[q][i] = prelim_composition_values[i][q];
cell->get_dof_indices (local_dof_indices);
@@ -845,12 +868,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);
Modified: branches/active_compositions/source/simulator/helper_functions.cc
===================================================================
--- branches/active_compositions/source/simulator/helper_functions.cc 2012-10-15 22:40:07 UTC (rev 1278)
+++ branches/active_compositions/source/simulator/helper_functions.cc 2012-10-16 18:05:48 UTC (rev 1279)
@@ -825,8 +825,17 @@
update_values | update_quadrature_points);
const FEValuesExtractors::Scalar pressure (dim);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int q=0;q<parameters.n_compositional_fields;++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ 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(),
@@ -840,14 +849,22 @@
fe_values.reinit (cell);
fe_values[pressure].get_function_values (this->solution,
pressure_values);
+ for(unsigned int i=0;i<parameters.n_compositional_fields;++i)
+ fe_values[compositional_fields[i]].get_function_values(this->solution,
+ composition_values[i]);
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;
@@ -889,8 +906,17 @@
update_quadrature_points );
const FEValuesExtractors::Scalar pressure (dim);
+ std::vector<FEValuesExtractors::Scalar> compositional_fields;
+ for (unsigned int q=0;q<parameters.n_compositional_fields;++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ 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(),
@@ -904,14 +930,22 @@
fe_values.reinit (cell);
fe_values[pressure].get_function_values (this->solution,
pressure_values);
+ for(unsigned int i=0;i<parameters.n_compositional_fields;++i)
+ fe_values[compositional_fields[i]].get_function_values(this->solution,
+ composition_values[i]);
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;
@@ -953,9 +987,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 q=0;q<parameters.n_compositional_fields;++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ 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(),
@@ -970,9 +1013,18 @@
pressure_values);
fe_values[temperature].get_function_values (this->solution,
temperature_values);
+ for(unsigned int i=0;i<parameters.n_compositional_fields;++i)
+ fe_values[compositional_fields[i]].get_function_values(this->solution,
+ composition_values[i]);
+
+ 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);
@@ -1010,9 +1062,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 q=0;q<parameters.n_compositional_fields;++q)
+ {
+ const FEValuesExtractors::Scalar temp(dim+2+q);
+ 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(),
@@ -1028,9 +1089,17 @@
pressure_values);
fe_values[temperature].get_function_values (this->solution,
temperature_values);
+ for(unsigned int i=0;i<parameters.n_compositional_fields;++i)
+ fe_values[compositional_fields[i]].get_function_values(this->solution,
+ composition_values[i]);
+ 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);
More information about the CIG-COMMITS
mailing list