From 6575547daf228bad0f3f16023fe5b41ed5e85a46 Mon Sep 17 00:00:00 2001 From: Ranpeng Li Date: Fri, 24 Apr 2026 17:57:12 +0200 Subject: [PATCH 1/2] specify entropy model's adiabatic profile via function --- doc/modules/changes/20260424_ranpengli | 4 + .../compute_entropy_profile.h | 22 ++++++ .../compute_entropy_profile.cc | 78 +++++++++++++++---- 3 files changed, 89 insertions(+), 15 deletions(-) create mode 100644 doc/modules/changes/20260424_ranpengli diff --git a/doc/modules/changes/20260424_ranpengli b/doc/modules/changes/20260424_ranpengli new file mode 100644 index 00000000000..c9679df8e93 --- /dev/null +++ b/doc/modules/changes/20260424_ranpengli @@ -0,0 +1,4 @@ +Added: An option to specify chemical compositions and their entropies +via functions when computing adiabatic profile of the entropy model. +
+(Ranpeng Li, 2026/04/24) diff --git a/include/aspect/adiabatic_conditions/compute_entropy_profile.h b/include/aspect/adiabatic_conditions/compute_entropy_profile.h index 277ab92b6f2..03a98e7a800 100644 --- a/include/aspect/adiabatic_conditions/compute_entropy_profile.h +++ b/include/aspect/adiabatic_conditions/compute_entropy_profile.h @@ -25,6 +25,7 @@ #include #include +#include namespace aspect { @@ -132,6 +133,27 @@ namespace aspect */ double delta_z; + /** + * An enum describing the different options to compute the reference + * profile for composition. + */ + enum CompositionProfile + { + surface_constant, + reference_function + }; + + /** + * Selected option to compute the reference profile for composition. + */ + CompositionProfile reference_composition; + + /** + * Function object that computes the reference composition profile + * if the reference_composition variable is set to function. + */ + std::unique_ptr> composition_function; + /** * A shared pointer to the initial composition object * that ensures that the current object can continue diff --git a/source/adiabatic_conditions/compute_entropy_profile.cc b/source/adiabatic_conditions/compute_entropy_profile.cc index d50c072e095..419233b8f87 100644 --- a/source/adiabatic_conditions/compute_entropy_profile.cc +++ b/source/adiabatic_conditions/compute_entropy_profile.cc @@ -96,21 +96,6 @@ namespace aspect // the entropy material model, and will therefore not be affected by this temperature. in.temperature[0] = this->get_adiabatic_surface_temperature(); - // Set all chemical composition to the initial composition, except the entropies, which - // are set to the surface entropy (since entropy is constant along an adiabat). - // Note, that if there a multiple entropy components they could have different entropies. - // However, since we are only interested in setting the - // equilibrated entropy, we do not need to compute the individual entropies for all components, - // and instead set all components to the equilibrated value. - // TODO : provide more ways to specify compositional fields like in compute_profile.cc - for (unsigned int c=0; cn_compositional_fields(); ++c) - { - if (this->introspection().get_composition_descriptions()[c].type == CompositionalFieldDescription::entropy) - in.composition[0][c] = surface_entropy; - else - in.composition[0][c] = initial_composition_manager->initial_composition(this->get_geometry_model().representative_point(0), c); - } - // Check whether gravity is pointing up / out or down / in. In the normal case it should // point down / in and therefore gravity should be positive, leading to increasing // adiabatic pressures and temperatures with depth. In some cases it will point up / out @@ -151,6 +136,32 @@ namespace aspect in.position[0] = representative_point; in.pressure[0] = pressures[i]; + + for (unsigned int c=0; cn_compositional_fields(); ++c) + { + // Set all chemical composition to the initial composition, except the entropies, which + // are set to the surface entropy (since entropy is constant along an adiabat). + // If there are multiple entropy components, which could have different entropies, + // we set one bulk entropy (which is the surface entropy) and calculate the equilibrated + // profile using the multi-component entropy method. + if (reference_composition == surface_constant) + if (this->introspection().get_composition_descriptions()[c].type == CompositionalFieldDescription::entropy) + in.composition[0][c] = surface_entropy; + else + in.composition[0][c] = initial_composition_manager->initial_composition(this->get_geometry_model().representative_point(0), c); + + // The chemical compositions and component's entropies can also be specified by a function. Each components + // will be equiliabrated using the multi-component entropy method. + else if (reference_composition == reference_function) + { + const Point<1> p(representative_point[1]); + for (unsigned int c=0; cn_compositional_fields(); ++c) + in.composition[0][c] = composition_function->value(p, c); + } + else + AssertThrow(false,ExcNotImplemented()); + } + this->get_material_model().evaluate(in, out); densities[i] = out.densities[0]; @@ -300,6 +311,14 @@ namespace aspect prm.declare_entry ("Surface entropy", "0", Patterns::Double(), "The surface entropy for the profile."); + + Functions::ParsedFunction<1>::declare_parameters (prm, 1); + prm.declare_entry ("Composition profile", "surface constant", + Patterns::Selection("surface constant|function"), + " Select the composition profile that is used to compute the" + " depth-dependent adiabatic profile. The default 'surface constant'" + " use the composition at the surface for all depth. You can instead " + " use 'function' to specify a depth-dependent composition variation"); } prm.leave_subsection(); } @@ -317,6 +336,35 @@ namespace aspect { n_points = prm.get_integer ("Number of points"); surface_entropy = prm.get_double ("Surface entropy"); + + const std::string composition_profile = prm.get("Composition profile"); + + if (composition_profile == "surface constant") + reference_composition = surface_constant; + else if (composition_profile == "function") + reference_composition = reference_function; + else + AssertThrow(false, ExcNotImplemented()); + + if ((this->n_compositional_fields() > 0) && (reference_composition == reference_function)) + { + composition_function + = std::make_unique>(this->n_compositional_fields()); + try + { + composition_function->parse_parameters (prm); + } + catch (...) + { + std::cerr << "ERROR: FunctionParser failed to parse\n" + << "\t'Adiabatic conditions model.compute entropy profile'\n" + << "with expression\n" + << "\t'" << prm.get("Function expression") << "'" + << "More information about the cause of the parse error \n" + << "is shown below.\n"; + throw; + } + } } prm.leave_subsection(); } From 6be2a0d0aab4c3fc09e3e74a698a3d59c1cec322 Mon Sep 17 00:00:00 2001 From: Ranpeng Li Date: Fri, 17 Jul 2026 16:10:35 +0200 Subject: [PATCH 2/2] changelog --- doc/modules/changes/20260717_ranpengli | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 doc/modules/changes/20260717_ranpengli diff --git a/doc/modules/changes/20260717_ranpengli b/doc/modules/changes/20260717_ranpengli new file mode 100644 index 00000000000..b96711de9f2 --- /dev/null +++ b/doc/modules/changes/20260717_ranpengli @@ -0,0 +1,4 @@ +Added: An option to specify chemical compositions and their entropies +via functions when computing adiabatic profile of the entropy model. +
+(Ranpeng Li, 2026/07/17)