Skip to content

Commit

Permalink
Extension to dynamic topography postprocessor to also calculate the t…
Browse files Browse the repository at this point in the history 
…opography at the lower surface; Adjust tests
  • Loading branch information
Jacqueline Austermann committed May 10, 2017
1 parent 86a9e02 commit 33c6c09
Show file tree
Hide file tree
Showing 19 changed files with 1,700 additions and 1,412 deletions.
4 changes: 4 additions & 0 deletions doc/modules/changes/20170508_austermann
View file
@@ -0,0 +1,4 @@
New: The dynamic topography post processor now also calculates the
topography at the bottom of the domain (and not only the upper surface).
<br>
(Jacky Austermann, 2017/05/08)
10 changes: 10 additions & 0 deletions include/aspect/postprocess/dynamic_topography.h
View file
Expand Up @@ -69,10 +69,20 @@ namespace aspect
*/
bool subtract_mean_dyn_topography;

/**
* A parameter that specifies whether the lower topography is outputted.
*/
bool output_lower_surface_topography;

/**
* A parameter allows users to set the density value outside the surface
*/
double density_above;

/**
* A parameter allows users to set the density value below the surface
*/
double density_below;
};
}
}
Expand Down
5 changes: 5 additions & 0 deletions include/aspect/postprocess/visualization/dynamic_topography.h
View file
Expand Up @@ -94,6 +94,11 @@ namespace aspect
* A parameter allows users to set the density value outside the surface
*/
double density_above;

/**
* A parameter allows users to set the density value below the lower surface
*/
double density_below;
};
}
}
Expand Down
218 changes: 177 additions & 41 deletions source/postprocess/dynamic_topography.cc
View file

Large diffs are not rendered by default.

117 changes: 90 additions & 27 deletions source/postprocess/visualization/dynamic_topography.cc
View file
Expand Up @@ -61,8 +61,10 @@ namespace aspect

std::vector<std::vector<double> > composition_values (this->n_compositional_fields(),std::vector<double> (quadrature_formula.size()));

double integrated_topography = 0;
double integrated_surface_area = 0;
double integrated_uppersurface_topography = 0;
double integrated_uppersurface_area = 0;
double integrated_lowersurface_topography = 0;
double integrated_lowersurface_area = 0;

// loop over all of the surface cells and if one less than h/3 away from
// the top surface, evaluate the stress at its center
Expand All @@ -75,21 +77,47 @@ namespace aspect
if (cell->is_locally_owned())
if (cell->at_boundary())
{
// see if the cell is at the *top* boundary, not just any boundary
unsigned int top_face_idx = numbers::invalid_unsigned_int;
{
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->at_boundary(f) && this->get_geometry_model().depth (cell->face(f)->center()) < cell->face(f)->minimum_vertex_distance()/3)
// see if the cell is at the *top* or *bottom* boundary, not just any boundary
unsigned int face_idx = numbers::invalid_unsigned_int;

// if the face is at the upper surface 'at_upper_surface' will be true, if
// it is at the lower surface 'at_upper_surface' will be false. The default
// is true and will be changed to false if it's at the lower boundary. If the
// cell is at neither boundary the loop will continue to the next cell.
bool at_upper_surface = true;
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
const double depth_face_center = this->get_geometry_model().depth (cell->face(f)->center());
const double upper_depth_cutoff = cell->face(f)->minimum_vertex_distance()/3.0;
const double lower_depth_cutoff = this->get_geometry_model().maximal_depth() - cell->face(f)->minimum_vertex_distance()/3.0;

// Check if cell is at upper and lower surface at the same time
if (depth_face_center < upper_depth_cutoff && depth_face_center > lower_depth_cutoff)
AssertThrow(false, ExcMessage("Your geometry model is so small that the upper and lower boundary of "
"the domain are bordered by the same cell. "
"Consider using a higher mesh resolution.") );

// Check if the face is at the top boundary
if (depth_face_center < upper_depth_cutoff)
{
top_face_idx = f;
face_idx = f;
break;
}
}
if (top_face_idx == numbers::invalid_unsigned_int)
// or at the bottom boundary
else if (depth_face_center > lower_depth_cutoff)
{
face_idx = f;
at_upper_surface = false;
break;
}
}

if (face_idx == numbers::invalid_unsigned_int)
{
(*return_value.second)(cell_index) = 0;
continue;
}

fe_values.reinit (cell);

// get the various components of the solution, then
Expand Down Expand Up @@ -123,7 +151,6 @@ 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;

Expand All @@ -145,7 +172,7 @@ namespace aspect
const double viscosity = out.viscosities[q];
const double density = out.densities[q];

const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q] - 1./3 * trace(in.strain_rate[q]) * unit_symmetric_tensor<dim>();
const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q] - 1./3.0 * trace(in.strain_rate[q]) * unit_symmetric_tensor<dim>();
const SymmetricTensor<2,dim> shear_stress = 2 * viscosity * strain_rate;

const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(location);
Expand All @@ -154,7 +181,11 @@ namespace aspect
// Subtract the dynamic pressure
const double dynamic_pressure = in.pressure[q] - this->get_adiabatic_conditions().pressure(location);
const double sigma_rr = gravity_direction * (shear_stress * gravity_direction) - dynamic_pressure;
const double dynamic_topography = - sigma_rr / (gravity_direction_binary*gravity.norm()) / (density - density_above);
// Choose the right density contrast depending on whether we're at the top or bottom boundary
const double density_outside = at_upper_surface ?
density_above :
density_below;
const double dynamic_topography = - sigma_rr / (gravity_direction_binary*gravity.norm()) / (density - density_outside);

// 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.
Expand All @@ -164,22 +195,32 @@ namespace aspect

const double dynamic_topography_cell_average = dynamic_topography_x_volume / volume;
// Compute the associated surface area to later compute the surfaces weighted integral
fe_face_values.reinit(cell, top_face_idx);
fe_face_values.reinit(cell, face_idx);
double surface = 0;
for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
{
surface += fe_face_values.JxW(q);
}

integrated_topography += dynamic_topography_cell_average*surface;
integrated_surface_area += surface;
if (at_upper_surface == true)
{
integrated_uppersurface_topography += dynamic_topography_cell_average*surface;
integrated_uppersurface_area += surface;
}
else
{
integrated_lowersurface_topography += dynamic_topography_cell_average*surface;
integrated_lowersurface_area += surface;
}

(*return_value.second)(cell_index) = dynamic_topography_cell_average;
}

// 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());
const double average_upper_topography = Utilities::MPI::sum (integrated_uppersurface_topography,this->get_mpi_communicator())
/ Utilities::MPI::sum (integrated_uppersurface_area,this->get_mpi_communicator());
const double average_lower_topography = Utilities::MPI::sum (integrated_lowersurface_topography,this->get_mpi_communicator())
/ Utilities::MPI::sum (integrated_lowersurface_area,this->get_mpi_communicator());

// subtract the average dynamic topography
cell_index = 0;
Expand All @@ -189,14 +230,22 @@ namespace aspect
if (cell->is_locally_owned()
&& (*return_value.second)(cell_index) != 0
&& cell->at_boundary())
{
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
if (cell->at_boundary(f) && this->get_geometry_model().depth (cell->face(f)->center()) < cell->face(f)->minimum_vertex_distance()/3)
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
const double depth_face_center = this->get_geometry_model().depth (cell->face(f)->center());
const double upper_depth_cutoff = cell->face(f)->minimum_vertex_distance()/3.0;
const double lower_depth_cutoff = this->get_geometry_model().maximal_depth() - cell->face(f)->minimum_vertex_distance()/3.0;
if (depth_face_center < upper_depth_cutoff)
{
(*return_value.second)(cell_index) -= average_topography;
(*return_value.second)(cell_index) -= average_upper_topography;
break;
}
}
else if (depth_face_center > lower_depth_cutoff)
{
(*return_value.second)(cell_index) -= average_lower_topography;
break;
}
}

return return_value;
}
Expand All @@ -215,9 +264,10 @@ namespace aspect
prm.declare_entry ("Subtract mean of dynamic topography", "false",
Patterns::Bool (),
"Option to remove the mean dynamic topography "
"in the outputted data file (not visualization). "
"in the visualization. The mean dynamic topography is calculated "
"and subtracted separately for the upper and lower boundary. "
"'true' subtracts the mean, 'false' leaves "
"the calculated dynamic topography as is. ");
"the calculated dynamic topography as is.");
prm.declare_entry ("Density above","0",
Patterns::Double (0),
"Dynamic topography is calculated as the excess or lack of mass that is supported by mantle flow. "
Expand All @@ -227,6 +277,15 @@ namespace aspect
"value of material that is displaced above the solid surface. By default this material is assumed to "
"be air, with a density of 0. "
"Units: $kg/m^3$.");
prm.declare_entry ("Density below","9900",
Patterns::Double (0),
"Dynamic topography is calculated as the excess or lack of mass that is supported by mantle flow. "
"This value depends on the density of material that is moved up or down, i.e. mantle above CMB, and the "
"density of the material that is displaced (generally outer core material). While the density of mantle rock "
"is part of the material model, this parameter `Density below' allows the user to specify the density "
"value of material that is displaced below the solid surface. By default this material is assumed to "
"be outer core material with a density of 9900. "
"Units: $kg/m^3$.");
}
prm.leave_subsection();
}
Expand All @@ -248,6 +307,7 @@ namespace aspect
{
subtract_mean_dyn_topography = prm.get_bool("Subtract mean of dynamic topography");
density_above = prm.get_double ("Density above");
density_below = prm.get_double ("Density below");
}
prm.leave_subsection();
}
Expand All @@ -270,7 +330,7 @@ namespace aspect
ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(DynamicTopography,
"dynamic topography",
"A visualization output object that generates output "
"for the dynamic topography. The approach to determine the "
"for the dynamic topography at the top and bottom of the model space. The approach to determine the "
"dynamic topography requires us to compute the stress tensor and "
"evaluate the component of it in the direction in which "
"gravity acts. In other words, we compute "
Expand All @@ -282,7 +342,10 @@ namespace aspect
"solve. From this, the dynamic "
"topography is computed using the formula "
"$h=\\frac{\\sigma_{rr}}{(\\mathbf g \\cdot \\mathbf n) \\rho}$ where $\\rho$ "
"is the density at the cell center. Note that this implementation takes "
"is the density at the cell center. For the bottom surface we chose the convection "
"that positive values are up (out) and negative values are in (down), analogous to "
"the deformation of the upper surface. "
"Note that this implementation takes "
"the direction of gravity into account, which means that reversing the flow "
"in backward advection calculations will not reverse the intantaneous topography "
"because the reverse flow will be divided by the reverse surface gravity."
Expand Down
2 changes: 1 addition & 1 deletion tests/dynamic_topography/dynamic_topography.00000
View file
@@ -1,4 +1,4 @@
# x y topography
# x y upper surface topography
112500 700000 -40.1439
337500 700000 -42.7954
562500 700000 -12.0079
Expand Down

0 comments on commit 33c6c09

Please sign in to comment.