diff --git a/include/aspect/mesh_deformation/fastscape.h b/include/aspect/mesh_deformation/fastscape.h index b06edad6cda..ee89efe0a34 100644 --- a/include/aspect/mesh_deformation/fastscape.h +++ b/include/aspect/mesh_deformation/fastscape.h @@ -453,6 +453,11 @@ namespace aspect */ double sediment_deposition_g; + /** + * Flag for allowing chemical compositions to have different erosional parameters for bedrock + */ + bool use_compositional_erosion_bedrock; + /** * Function of bedrock river incision rate (kf) for the stream power law. * Represents the parameter `kf` in the FastScape landscape evolution equation. @@ -475,7 +480,7 @@ namespace aspect * otherwise, the units are ${m^(1-2drainage_area_exponent)/s}$. Then a time scale factor will be applied to * convert it into ${m^(1-2drainage_area_exponent)/yr}$ for Fastscape. */ - double constant_bedrock_river_incision_rate; + std::vector constant_bedrock_river_incision_rate; /** * Sediment river incision rate for the stream power law. @@ -509,7 +514,7 @@ namespace aspect * convert it into ${m^2/yr}$ for Fastscape. * This function is only used only if use_kf_distribution_function is false. */ - double constant_bedrock_transport_coefficient; + std::vector constant_bedrock_transport_coefficient; /** * Sediment transport coefficient for hillslope diffusion. diff --git a/source/mesh_deformation/fastscape.cc b/source/mesh_deformation/fastscape.cc index 58d2f6838e6..80ec38fbc16 100644 --- a/source/mesh_deformation/fastscape.cc +++ b/source/mesh_deformation/fastscape.cc @@ -526,7 +526,7 @@ namespace aspect // (by default -1). If a different rate is set for sediment and // bedrock, the sediment rate is only used for sediment layers // at least 1 m thick. - if (sediment_river_incision_rate < 0. || sediment_thickness <= 1.) + if (sediment_river_incision_rate < 0. || sediment_thickness <= 1. || std::isnan(sediment_thickness)) combined_kf[i] = bedrock_river_incision_rate_array[i]; else if (sediment_river_incision_rate >= 0. && sediment_thickness > 1.) combined_kf[i] = sediment_river_incision_rate; @@ -540,7 +540,7 @@ namespace aspect // (by default -1). If a different rate is set for sediment and bedrock, // the sediment coefficient is only used for sediment layers // at least 1 m thick. - if (sediment_transport_coefficient < 0 || sediment_thickness <= 1.) + if (sediment_transport_coefficient < 0 || sediment_thickness <= 1. || std::isnan(sediment_thickness)) combined_kd[i] = bedrock_transport_coefficient_array[i]; else if (sediment_transport_coefficient >= 0. && sediment_thickness > 1.) combined_kd[i] = sediment_transport_coefficient; @@ -685,8 +685,12 @@ namespace aspect FastScape::get_aspect_values() const { + std::vector compositional_field_names = this->introspection().get_composition_names(); + const std::vector chemical_composition_idx = this->introspection().chemical_composition_field_indices(); + const unsigned int n_chemical_compositional_fields = chemical_composition_idx.size(); + const types::boundary_id relevant_boundary = this->get_geometry_model().translate_symbolic_boundary_name_to_id ("top"); - std::vector> local_aspect_values(dim+2, std::vector()); + std::vector> local_aspect_values(dim+4, std::vector()); // Get a quadrature rule that exists only on the corners, and increase the refinement if specified. const QIterated face_corners (QTrapezoid<1>(), @@ -706,10 +710,23 @@ namespace aspect if ( cell->face(face_no)->boundary_id() != relevant_boundary) continue; + std::vector> composition_values_array(this->n_compositional_fields(), std::vector(face_corners.size())); std::vector> vel(face_corners.size()); fe_face_values.reinit(cell, face_no); fe_face_values[this->introspection().extractors.velocities].get_function_values(this->get_solution(), vel); + if (use_compositional_erosion_bedrock) + { + // Get compositional erosional parameters from compositional values, for chemical compositions + background mantle + for (unsigned int c=0; cintrospection().extractors.compositional_fields[field]].get_function_values(this->get_solution(), composition_values_array[field]); + Assert(composition_values_array[field].size() > 0, + ExcMessage("Composition_values_array[c] is empty")); + } + } + for (unsigned int corner = 0; corner < face_corners.size(); ++corner) { const Point vertex = fe_face_values.quadrature_point(corner); @@ -726,6 +743,22 @@ namespace aspect if (std::abs(indx - std::round(indx)) >= node_tolerance) continue; + double bedrock_river_incision_rate_at_point = numbers::signaling_nan(); + double bedrock_transport_coefficient_at_point = numbers::signaling_nan(); + if (use_compositional_erosion_bedrock) + { + std::vector composition_values(n_chemical_compositional_fields); + double volume_fraction_sum = 0; + for (unsigned int c=0; c < n_chemical_compositional_fields; ++c) + { + composition_values[c] = composition_values_array[chemical_composition_idx[c]][corner]; + volume_fraction_sum = volume_fraction_sum + composition_values_array[chemical_composition_idx[c]][corner]; + } + // Add background material fraction at the beginning + composition_values.insert(composition_values.begin(), std::max(0.0, 1.0 - volume_fraction_sum)); + bedrock_river_incision_rate_at_point = MaterialModel::MaterialUtilities::average_value (composition_values, constant_bedrock_river_incision_rate, MaterialModel::MaterialUtilities::arithmetic); + bedrock_transport_coefficient_at_point = MaterialModel::MaterialUtilities::average_value (composition_values, constant_bedrock_transport_coefficient, MaterialModel::MaterialUtilities::arithmetic); + } // If we're in 2D, we want to take the values and apply them to every row of X points. if (dim == 2) @@ -752,6 +785,9 @@ namespace aspect // Always convert to m/yr for FastScape local_aspect_values[2+d].push_back(vel[corner][d]*year_in_seconds); } + + local_aspect_values[dim+2].push_back(bedrock_river_incision_rate_at_point); + local_aspect_values[dim+3].push_back(bedrock_transport_coefficient_at_point); } } // 3D case @@ -776,6 +812,9 @@ namespace aspect { local_aspect_values[2+d].push_back(vel[corner][d]*year_in_seconds); } + + local_aspect_values[dim+2].push_back(bedrock_river_incision_rate_at_point); + local_aspect_values[dim+3].push_back(bedrock_transport_coefficient_at_point); } } } @@ -793,6 +832,8 @@ namespace aspect std::vector &velocity_z, std::vector> &local_aspect_values) const { + const double time_scaling_factor = (this->convert_output_to_years() ? 1.0 : year_in_seconds); + for (unsigned int i=0; iget_mpi_communicator()); ++p) @@ -843,14 +890,23 @@ namespace aspect velocity_y[index] = 0; else velocity_y[index] = local_aspect_values[3][i]; + + if (use_compositional_erosion_bedrock) + { + bedrock_river_incision_rate_array[index] = time_scaling_factor * local_aspect_values[dim+2][i]; + bedrock_transport_coefficient_array[index] = time_scaling_factor * local_aspect_values[dim+3][i]; + } + } } + bool fastscape_mesh_filled = true; + // Initialize the bedrock river incision rate and transport coefficient, // and check that there are no empty mesh points due to // an improperly set maximum_surface_refinement_level, additional_refinement_levels, // and surface_refinement_difference - bool fastscape_mesh_filled = true; + const unsigned int fastscape_array_size = fastscape_nx*fastscape_ny; for (unsigned int i=0; iconvert_output_to_years() ? 1.0 : year_in_seconds); - // Set bedrock transport coefficient kd either from a function or a constant. - bedrock_transport_coefficient_array[i] = - (use_kd_distribution_function - ? - time_scaling_factor * kd_distribution_function.value(Point<2>(x, y)) - : - time_scaling_factor * constant_bedrock_transport_coefficient); - - // Set bedrock river incision rate kf either from a function or a constant. - bedrock_river_incision_rate_array[i] = - (use_kf_distribution_function) - ? - time_scaling_factor * kf_distribution_function.value(Point<2>(x, y)) - : - time_scaling_factor * constant_bedrock_river_incision_rate; - - - // If this is a boundary node that is a ghost node then ignore that it - // has not filled yet as the ghost nodes haven't been set. - if (elevation[i] == std::numeric_limits::max() && !is_ghost_node(i,false)) - fastscape_mesh_filled = false; + if (!use_kf_distribution_function) + { + bedrock_river_incision_rate_array[i] = time_scaling_factor * constant_bedrock_river_incision_rate[0]; + } + else + bedrock_river_incision_rate_array[i] = time_scaling_factor * kf_distribution_function.value(Point<2>(x, y)); + if (!use_kd_distribution_function) + { + bedrock_transport_coefficient_array[i] = time_scaling_factor * constant_bedrock_transport_coefficient[0]; + } + else + bedrock_transport_coefficient_array[i] = time_scaling_factor * kd_distribution_function.value(Point<2>(x, y)); + + if (elevation[i] == std::numeric_limits::max() && !is_ghost_node(i,false) || std::isnan(elevation[i])) + { + fastscape_mesh_filled = false; + } } + // If this is a boundary node that is a ghost node then ignore that it + // has not filled yet as the ghost nodes haven't been set. fastscape_mesh_filled = Utilities::MPI::broadcast(this->get_mpi_communicator(), fastscape_mesh_filled, 0); AssertThrow (fastscape_mesh_filled == true, ExcMessage("The FastScape mesh is missing data. A likely cause for this is that the " @@ -1837,6 +1884,10 @@ namespace aspect Patterns::Double(), "Deposition coefficient for sediment. A value smaller than 0 sets this to the same as the bedrock deposition coefficient."); // Define Bedrock river incision rate (Kf) as a constant value or a time-dependent user-defined function + prm.declare_entry("Allow compositional erosion for bedrock", "false", + Patterns::Bool (), + "Whether to allow chemical compositions to have different erosional parameters for bedrock, " + "including river incision rate and diffusivity."); prm.declare_entry("Use kf distribution function", "false", Patterns::Bool(), "Whether to define bedrock river incision rate using a distribution function. " @@ -1844,8 +1895,9 @@ namespace aspect "the parameter ``Bedrock river incision rate''. Units: ${m^(1-2drainage_area_exponent)/yr}$ " "if ``Use years instead of seconds'' is true; otherwise, the units are ${m^(1-2drainage_area_exponent)/s}$."); prm.declare_entry("Bedrock river incision rate", "1e-5", - Patterns::Double(), - "River incision rate for bedrock in the Stream Power Law. " + Patterns::List(Patterns::Double(0.)), + "A list of river incision rate for bedrock in the Stream Power Law, for background mantle and chemical " + "composition fields, for a total of N+1 values, where N is the number of chemical composition fields. " "Units: ${m^(1-2drainage_area_exponent)/yr}$ if ``Use years instead of seconds'' is true; " "otherwise, the units are ${m^(1-2drainage_area_exponent)/s}$."); prm.enter_subsection ("kf distribution function"); @@ -1867,9 +1919,10 @@ namespace aspect "``Bedrock diffusivity''. Units: ${m^2/yr}$ if ``Use years instead of seconds'' " "is true; otherwise, the units are ${m^2/s}$."); prm.declare_entry("Bedrock diffusivity", "1e-2", - Patterns::Double(), - "Transport coefficient (diffusivity) for bedrock. Units: ${m^2/yr}$ if ``Use years instead of seconds'' " - "is true; otherwise, the units are ${m^2/s}$."); + Patterns::List(Patterns::Double(0.)), + "A list of transport coefficient (diffusivity) for bedrock, for background mantle and chemical " + "composition fields, for a total of N+1 values, where N is the number of chemical composition fields. " + " Units: ${m^2/yr}$ if ``Use years instead of seconds'' is true; otherwise, the units are ${m^2/s}$."); prm.enter_subsection ("kd distribution function"); { Functions::ParsedFunction<2>::declare_parameters(prm, 2); @@ -1962,7 +2015,7 @@ namespace aspect prm.declare_entry("Depth averaging thickness", "1e2", Patterns::Double(), "Depth averaging for the sand-silt equation. Units: ${m}$"); - prm.declare_entry("Sand transport coefficient", "5e2", + prm.declare_entry("Sand transport coefficient", "2.5e2", Patterns::Double(), "Transport coefficient (diffusivity) for sand. Units: ${m^2/yr}$"); prm.declare_entry("Silt transport coefficient", "2.5e2", @@ -2070,8 +2123,17 @@ namespace aspect // with a year in seconds. The bedrock values are scaled when filling the FastScape // arrays. const double time_scaling_factor = (this->convert_output_to_years() ? 1.0 : year_in_seconds); + sediment_river_incision_rate = time_scaling_factor * prm.get_double("Sediment river incision rate"); + use_compositional_erosion_bedrock = prm.get_bool("Allow compositional erosion for bedrock"); // kf + // Make options file for parsing maps to double arrays for bedrock river incision rate and bedrock transport coefficient + std::vector chemical_field_names = this->introspection().chemical_composition_field_names(); + // Establish that a background field is required here + chemical_field_names.insert(chemical_field_names.begin(),"background"); + Utilities::MapParsing::Options options(chemical_field_names, ""); + options.list_of_allowed_keys = chemical_field_names; + options.property_name = "Bedrock river incision rate"; use_kf_distribution_function = prm.get_bool("Use kf distribution function"); if (use_kf_distribution_function) { @@ -2083,11 +2145,14 @@ namespace aspect // If using the function description for the kf, no parts of the code // base should use the constant_bedrock_river_incision_rate variable. Poison it to // make sure it really isn't used anywhere: - constant_bedrock_river_incision_rate = numbers::signaling_nan(); + constant_bedrock_river_incision_rate.assign( + constant_bedrock_river_incision_rate.size(), + numbers::signaling_nan() + ); } else { - constant_bedrock_river_incision_rate = prm.get_double("Bedrock river incision rate"); + constant_bedrock_river_incision_rate = Utilities::MapParsing::parse_map_to_double_array(prm.get(options.property_name), options); } sediment_transport_coefficient = time_scaling_factor * prm.get_double("Sediment diffusivity"); // kd @@ -2102,12 +2167,19 @@ namespace aspect // If using the function description for the kd, no parts of the code // base should use the constant_bedrock_river_incision_rate variable. Poison it to // make sure it really isn't used anywhere: - constant_bedrock_transport_coefficient = numbers::signaling_nan(); + constant_bedrock_transport_coefficient.assign( + constant_bedrock_transport_coefficient.size(), + numbers::signaling_nan() + ); } else { - constant_bedrock_transport_coefficient = prm.get_double("Bedrock diffusivity"); + options.property_name = "Bedrock diffusivity"; + constant_bedrock_transport_coefficient = Utilities::MapParsing::parse_map_to_double_array(prm.get(options.property_name), options); } + AssertThrow(!(use_compositional_erosion_bedrock && + (use_kf_distribution_function || use_kd_distribution_function)), + ExcMessage("Compositional erosion for bedrock cannot be used together with kf or kd distribution functions.")); bedrock_deposition_g = prm.get_double("Bedrock deposition coefficient"); sediment_deposition_g = prm.get_double("Sediment deposition coefficient"); slope_exponent_p = prm.get_double("Multi-direction slope exponent"); @@ -2168,7 +2240,6 @@ namespace aspect sand_silt_averaging_depth = prm.get_double("Depth averaging thickness"); sand_transport_coefficient = prm.get_double("Sand transport coefficient"); silt_transport_coefficient = prm.get_double("Silt transport coefficient"); - if (!this->convert_output_to_years()) { sand_transport_coefficient *= year_in_seconds;