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
  • Loading branch information
Jacqueline Austermann committed May 9, 2017
1 parent 4cf330c commit a519f48
Show file tree
Hide file tree
Showing 5 changed files with 284 additions and 71 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: 175 additions & 43 deletions source/postprocess/dynamic_topography.cc
View file
Expand Up @@ -54,14 +54,18 @@ namespace aspect

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

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

// have a stream into which we write the data. the text stream is then
// later sent to processor 0
std::ostringstream output;

double integrated_topography = 0;
double integrated_surface_area = 0;
std::ostringstream output_upper;
std::ostringstream output_lower;

std::vector<std::pair<Point<dim>,double> > stored_values;
std::vector<std::pair<Point<dim>,double> > stored_values_upper;
std::vector<std::pair<Point<dim>,double> > stored_values_lower;

// 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 @@ -73,19 +77,43 @@ 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
bool at_upper_surface;
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 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)
{
face_idx = f;
at_upper_surface = true;
break;
}
// or at the bottom boundary
else if (depth_face_center > lower_depth_cutoff && output_lower_surface_topography == true)
{
top_face_idx = f;
face_idx = f;
at_upper_surface = false;
break;
}
}

if (face_idx == numbers::invalid_unsigned_int)
continue;

if (top_face_idx == numbers::invalid_unsigned_int)
continue;
}
fe_values.reinit (cell);

// get the various components of the solution, then
Expand Down Expand Up @@ -142,7 +170,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 @@ -151,7 +179,12 @@ 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 @@ -161,39 +194,49 @@ 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);
}

const Point<dim> midpoint_at_surface = cell->face(top_face_idx)->center();
const Point<dim> midpoint_at_surface = cell->face(face_idx)->center();

integrated_topography += dynamic_topography_cell_average*surface;
integrated_surface_area += surface;

stored_values.push_back (std::make_pair(midpoint_at_surface, dynamic_topography_cell_average));
if (at_upper_surface == true)
{
integrated_uppersurface_topography += dynamic_topography_cell_average*surface;
integrated_uppersurface_area += surface;
stored_values_upper.push_back (std::make_pair(midpoint_at_surface, dynamic_topography_cell_average));
}
else
{
integrated_lowersurface_topography += dynamic_topography_cell_average*surface;
integrated_lowersurface_area += surface;
stored_values_lower.push_back (std::make_pair(midpoint_at_surface, 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());

// Write the solution to an output file
// if (DT_mean_switch == true) subtract the average dynamic topography,
// otherwise leave as is
for (unsigned int i=0; i<stored_values.size(); ++i)
for (unsigned int i=0; i<stored_values_upper.size(); ++i)
{
output << stored_values[i].first
<< ' '
<< stored_values[i].second -
(subtract_mean_dyn_topography
?
average_topography
:
0.)
<< std::endl;
output_upper << stored_values_upper[i].first
<< ' '
<< stored_values_upper[i].second -
(subtract_mean_dyn_topography
?
average_upper_topography
:
0.)
<< std::endl;
}


Expand All @@ -203,7 +246,7 @@ namespace aspect
if (this->get_parameters().run_postprocessors_on_nonlinear_iterations)
filename.append("." + Utilities::int_to_string (this->get_nonlinear_iteration(), 4));

const unsigned int max_data_length = Utilities::MPI::max (output.str().size()+1,
const unsigned int max_data_length = Utilities::MPI::max (output_upper.str().size()+1,
this->get_mpi_communicator());
const unsigned int mpi_tag = 123;

Expand All @@ -215,10 +258,10 @@ namespace aspect

file << "# "
<< ((dim==2)? "x y" : "x y z")
<< " topography" << std::endl;
<< " upper surface topography" << std::endl;

// first write out the data we have created locally
file << output.str();
file << output_upper.str();

std::string tmp;
tmp.resize (max_data_length, '\0');
Expand Down Expand Up @@ -247,10 +290,81 @@ namespace aspect
// on other processors, send the data to processor zero. include the \0
// character at the end of the string
{
MPI_Send (&output.str()[0], output.str().size()+1, MPI_CHAR, 0, mpi_tag,
MPI_Send (&output_upper.str()[0], output_upper.str().size()+1, MPI_CHAR, 0, mpi_tag,
this->get_mpi_communicator());
}


if (output_lower_surface_topography == true)
{
for (unsigned int i=0; i<stored_values_lower.size(); ++i)
{
output_lower << stored_values_lower[i].first
<< ' '
<< stored_values_lower[i].second -
(subtract_mean_dyn_topography
?
average_lower_topography
:
0.)
<< std::endl;
}

std::string filename = this->get_output_directory() +
"bottom_dynamic_topography." +
Utilities::int_to_string(this->get_timestep_number(), 5);
if (this->get_parameters().run_postprocessors_on_nonlinear_iterations)
filename.append("." + Utilities::int_to_string (this->get_nonlinear_iteration(), 4));

const unsigned int max_data_length = Utilities::MPI::max (output_lower.str().size()+1,
this->get_mpi_communicator());
const unsigned int mpi_tag = 123;

// on processor 0, collect all of the data the individual processors send
// and concatenate them into one file
if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0)
{
std::ofstream file (filename.c_str());

file << "# "
<< ((dim==2)? "x y" : "x y z")
<< " lower surface topography" << std::endl;

// first write out the data we have created locally
file << output_lower.str();

std::string tmp;
tmp.resize (max_data_length, '\0');

// then loop through all of the other processors and collect
// data, then write it to the file
for (unsigned int p=1; p<Utilities::MPI::n_mpi_processes(this->get_mpi_communicator()); ++p)
{
MPI_Status status;
// get the data. note that MPI says that an MPI_Recv may receive
// less data than the length specified here. since we have already
// determined the maximal message length, we use this feature here
// rather than trying to find out the exact message length with
// a call to MPI_Probe.
MPI_Recv (&tmp[0], max_data_length, MPI_CHAR, p, mpi_tag,
this->get_mpi_communicator(), &status);

// output the string. note that 'tmp' has length max_data_length,
// but we only wrote a certain piece of it in the MPI_Recv, ended
// by a \0 character. write only this part by outputting it as a
// C string object, rather than as a std::string
file << tmp.c_str();
}
}
else
// on other processors, send the data to processor zero. include the \0
// character at the end of the string
{
MPI_Send (&output_lower.str()[0], output_lower.str().size()+1, MPI_CHAR, 0, mpi_tag,
this->get_mpi_communicator());
}
}

return std::pair<std::string,std::string>("Writing dynamic topography:",
filename);
}
Expand All @@ -267,17 +381,31 @@ 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("Output lower surface topography","false",
Patterns::Bool (),
"Option to output the lower surface topography "
"in a separate file lower_dynamic_topography.");
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. "
"This value depends on the density of material that is moved up or down, i.e. crustal rock, and the "
"density of the material that is displaced (generally water or air). While the density of crustal rock "
"is part of the material model, this parameter `Density above' allows the user to specify the density "
"value of material that is displaced above the solid surface. By default this material is assumed to "
"be air, with a density of 0. "
"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 @@ -294,7 +422,9 @@ namespace aspect
prm.enter_subsection("Dynamic Topography");
{
subtract_mean_dyn_topography = prm.get_bool("Subtract mean of dynamic topography");
output_lower_surface_topography = prm.get_bool("Output lower surface topography");
density_above = prm.get_double ("Density above");
density_below = prm.get_double ("Density below");
}
prm.leave_subsection();
}
Expand All @@ -313,7 +443,7 @@ namespace aspect
ASPECT_REGISTER_POSTPROCESSOR(DynamicTopography,
"dynamic topography",
"A postprocessor that computes a measure of dynamic topography "
"based on the stress at the surface. The data is written into text "
"based on the stress at the surface and bottom. The data is written into text "
"files named 'dynamic\\_topography.NNNNN' in the output directory, "
"where NNNNN is the number of the time step."
"\n\n"
Expand All @@ -329,7 +459,9 @@ 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 instantaneous topography "
"because the reverse flow will be divided by the reverse surface gravity. "
Expand Down

0 comments on commit a519f48

Please sign in to comment.