[cig-commits] [commit] master: Change averaging in dynamic topography postprocessor from volume weighted to surface area weighted (14d2520)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Sun Jan 18 17:48:53 PST 2015
Repository : https://github.com/geodynamics/aspect
On branch : master
Link : https://github.com/geodynamics/aspect/compare/eef6299b61e1952be22819a80f6a6e4ee344b677...09cd1e0322ded4719e699729a269b602229b86b3
>---------------------------------------------------------------
commit 14d2520ee88e7521c7a84a734faa78fc48291b07
Author: Jacky Austermann <jaustermann at fas.harvard.edu>
Date: Wed Jan 14 18:59:48 2015 -0500
Change averaging in dynamic topography postprocessor from volume weighted to surface area weighted
>---------------------------------------------------------------
14d2520ee88e7521c7a84a734faa78fc48291b07
source/postprocess/dynamic_topography.cc | 39 ++++++++++++++++----
.../visualization/dynamic_topography.cc | 42 +++++++++++++++++-----
2 files changed, 67 insertions(+), 14 deletions(-)
diff --git a/source/postprocess/dynamic_topography.cc b/source/postprocess/dynamic_topography.cc
index 4ff4b34..f88260c 100644
--- a/source/postprocess/dynamic_topography.cc
+++ b/source/postprocess/dynamic_topography.cc
@@ -34,6 +34,7 @@ namespace aspect
DynamicTopography<dim>::execute (TableHandler &statistics)
{
const QMidpoint<dim> quadrature_formula;
+ const QMidpoint<dim-1> quadrature_formula_face;
FEValues<dim> fe_values (this->get_mapping(),
this->get_fe(),
@@ -41,6 +42,10 @@ namespace aspect
update_values |
update_gradients |
update_q_points);
+ FEFaceValues<dim> fe_face_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature_formula_face,
+ update_JxW_values);
typename MaterialModel::Interface<dim>::MaterialModelInputs in(fe_values.n_quadrature_points, this->n_compositional_fields());
typename MaterialModel::Interface<dim>::MaterialModelOutputs out(fe_values.n_quadrature_points, this->n_compositional_fields());
@@ -110,9 +115,14 @@ namespace aspect
// for each of the quadrature points, evaluate the
// stress and compute the component in direction of the
// gravity vector
+
+ double dynamic_topography_x_volume = 0;
+ double volume = 0;
+ Point<dim> location = fe_values.quadrature_point(0);
+
for (unsigned int q=0; q<quadrature_formula.size(); ++q)
{
- const Point<dim> location = fe_values.quadrature_point(q);
+ location = fe_values.quadrature_point(q);
const double viscosity = out.viscosities[q];
const double density = out.densities[q];
@@ -127,13 +137,30 @@ namespace aspect
const double sigma_rr = gravity_direction * (shear_stress * gravity_direction) - dynamic_pressure;
const double dynamic_topography = - sigma_rr / gravity.norm() / density;
- integrated_topography += dynamic_topography * fe_values.JxW(q);
- integrated_surface_area += fe_values.JxW(q);
-
- stored_values.push_back (std::make_pair(location, dynamic_topography));
- }
+ // JxW provides the volume quadrature weights. This is a general formulation
+ // necessary for when a quadrature formula is used that has more than one point.
+ dynamic_topography_x_volume += dynamic_topography * fe_values.JxW(q);
+ volume += fe_values.JxW(q);
+ }
+
+ const double dynamic_topography_total = dynamic_topography_x_volume / volume;
+ // Compute the associated surface area to later compute the surfaces weighted integral
+ double surface = 0;
+ for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->at_boundary(f))
+ if (this->get_geometry_model().depth (cell->face(f)->center()) < cell->face(f)->minimum_vertex_distance()/3)
+ {
+ fe_face_values.reinit(cell,f);
+ surface = fe_face_values.JxW(0);
+ }
+
+ integrated_topography += dynamic_topography_total*surface;
+ integrated_surface_area += surface;
+
+ stored_values.push_back (std::make_pair(location, dynamic_topography_total));
}
+ // Calculate surface weighted average dynamic topography
const double average_topography = Utilities::MPI::sum (integrated_topography,this->get_mpi_communicator()) / Utilities::MPI::sum (integrated_surface_area,this->get_mpi_communicator());
diff --git a/source/postprocess/visualization/dynamic_topography.cc b/source/postprocess/visualization/dynamic_topography.cc
index 0aff358..05d86e7 100644
--- a/source/postprocess/visualization/dynamic_topography.cc
+++ b/source/postprocess/visualization/dynamic_topography.cc
@@ -42,7 +42,7 @@ namespace aspect
// evaluate a single point per cell
const QMidpoint<dim> quadrature_formula;
- const unsigned int n_q_points = quadrature_formula.size();
+ const QMidpoint<dim-1> quadrature_formula_face;
FEValues<dim> fe_values (this->get_mapping(),
this->get_fe(),
@@ -51,6 +51,11 @@ namespace aspect
update_gradients |
update_quadrature_points );
+ FEFaceValues<dim> fe_face_values (this->get_mapping(),
+ this->get_fe(),
+ quadrature_formula_face,
+ update_JxW_values);
+
typename MaterialModel::Interface<dim>::MaterialModelInputs in(fe_values.n_quadrature_points, this->n_compositional_fields());
typename MaterialModel::Interface<dim>::MaterialModelOutputs out(fe_values.n_quadrature_points, this->n_compositional_fields());
@@ -117,6 +122,10 @@ namespace aspect
// for each of the quadrature points, evaluate the
// stress and compute the component in direction of the
// gravity vector
+
+ double dynamic_topography_x_volume = 0;
+ double volume = 0;
+
for (unsigned int q=0; q<quadrature_formula.size(); ++q)
{
const Point<dim> location = fe_values.quadrature_point(q);
@@ -134,13 +143,30 @@ namespace aspect
const double sigma_rr = gravity_direction * (shear_stress * gravity_direction) - dynamic_pressure;
const double dynamic_topography = - sigma_rr / gravity.norm() / density;
- integrated_topography += dynamic_topography * fe_values.JxW(q);
- integrated_surface_area += fe_values.JxW(q);
-
- (*return_value.second)(cell_index) = dynamic_topography;
- }
- }
-
+ // JxW provides the volume quadrature weights. This is a general formulation
+ // necessary for when a quadrature formula is used that has more than one point.
+ dynamic_topography_x_volume += dynamic_topography * fe_values.JxW(q);
+ volume += fe_values.JxW(q);
+ }
+
+ const double dynamic_topography_total = dynamic_topography_x_volume / volume;
+ // Compute the associated surface area to later compute the surfaces weighted integral
+ double surface = 0;
+ for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
+ if (cell->at_boundary(f))
+ if (this->get_geometry_model().depth (cell->face(f)->center()) < cell->face(f)->minimum_vertex_distance()/3)
+ {
+ fe_face_values.reinit(cell,f);
+ surface = fe_face_values.JxW(0);
+ }
+
+ integrated_topography += dynamic_topography_total*surface;
+ integrated_surface_area += surface;
+
+ (*return_value.second)(cell_index) = dynamic_topography_total;
+ }
+
+ // Calculate surface weighted average dynamic topography
const double average_topography = Utilities::MPI::sum (integrated_topography,this->get_mpi_communicator()) / Utilities::MPI::sum (integrated_surface_area,this->get_mpi_communicator());
// if (DT_mean_switch == true) subtract the average dynamic topography,
More information about the CIG-COMMITS
mailing list