From a519f489b69972048720405d9378cbc9e95cff4e Mon Sep 17 00:00:00 2001 From: Jacqueline Austermann Date: Mon, 8 May 2017 14:40:05 -0400 Subject: [PATCH] Extension to dynamic topography postprocessor to also calculate the topography at the lower surface --- doc/modules/changes/20170508_austermann | 4 + .../aspect/postprocess/dynamic_topography.h | 10 + .../visualization/dynamic_topography.h | 5 + source/postprocess/dynamic_topography.cc | 218 ++++++++++++++---- .../visualization/dynamic_topography.cc | 118 +++++++--- 5 files changed, 284 insertions(+), 71 deletions(-) create mode 100644 doc/modules/changes/20170508_austermann diff --git a/doc/modules/changes/20170508_austermann b/doc/modules/changes/20170508_austermann new file mode 100644 index 00000000000..35d6157cb6f --- /dev/null +++ b/doc/modules/changes/20170508_austermann @@ -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). +
+(Jacky Austermann, 2017/05/08) diff --git a/include/aspect/postprocess/dynamic_topography.h b/include/aspect/postprocess/dynamic_topography.h index bdbca9627ec..d72d25a97b5 100644 --- a/include/aspect/postprocess/dynamic_topography.h +++ b/include/aspect/postprocess/dynamic_topography.h @@ -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; }; } } diff --git a/include/aspect/postprocess/visualization/dynamic_topography.h b/include/aspect/postprocess/visualization/dynamic_topography.h index e1fd38f9a1c..8c062cc993d 100644 --- a/include/aspect/postprocess/visualization/dynamic_topography.h +++ b/include/aspect/postprocess/visualization/dynamic_topography.h @@ -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; }; } } diff --git a/source/postprocess/dynamic_topography.cc b/source/postprocess/dynamic_topography.cc index 22ae25c96a3..380b1676320 100644 --- a/source/postprocess/dynamic_topography.cc +++ b/source/postprocess/dynamic_topography.cc @@ -54,14 +54,18 @@ namespace aspect std::vector > composition_values (this->n_compositional_fields(),std::vector (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,double> > stored_values; + std::vector,double> > stored_values_upper; + std::vector,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 @@ -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::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::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 @@ -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(); + const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q] - 1./3.0 * trace(in.strain_rate[q]) * unit_symmetric_tensor(); const SymmetricTensor<2,dim> shear_stress = 2 * viscosity * strain_rate; const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(location); @@ -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. @@ -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 midpoint_at_surface = cell->face(top_face_idx)->center(); + const Point 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; iget_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; @@ -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'); @@ -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; iget_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; pget_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("Writing dynamic topography:", filename); } @@ -267,9 +381,14 @@ 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. " @@ -277,7 +396,16 @@ namespace aspect "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(); @@ -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(); } @@ -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" @@ -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. " diff --git a/source/postprocess/visualization/dynamic_topography.cc b/source/postprocess/visualization/dynamic_topography.cc index a5d3c8388ee..b185bb8ae25 100644 --- a/source/postprocess/visualization/dynamic_topography.cc +++ b/source/postprocess/visualization/dynamic_topography.cc @@ -61,8 +61,10 @@ namespace aspect std::vector > composition_values (this->n_compositional_fields(),std::vector (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 @@ -75,21 +77,46 @@ 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::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::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) { - top_face_idx = f; + face_idx = f; + at_upper_surface = true; 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 @@ -123,7 +150,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; @@ -145,7 +171,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(); + const SymmetricTensor<2,dim> strain_rate = in.strain_rate[q] - 1./3.0 * trace(in.strain_rate[q]) * unit_symmetric_tensor(); const SymmetricTensor<2,dim> shear_stress = 2 * viscosity * strain_rate; const Tensor<1,dim> gravity = this->get_gravity_model().gravity_vector(location); @@ -154,7 +180,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. @@ -164,22 +194,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; qget_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; @@ -189,14 +229,22 @@ namespace aspect if (cell->is_locally_owned() && (*return_value.second)(cell_index) != 0 && cell->at_boundary()) - { - for (unsigned int f=0; f::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::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; } @@ -215,9 +263,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. " @@ -225,7 +274,16 @@ namespace aspect "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(); @@ -248,6 +306,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(); } @@ -270,7 +329,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 " @@ -282,7 +341,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."