# Problems with initial composition model using reactive fluid transport

**URL:** <https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329>\
**Category:** ASPECT\
**Created:** [October 16, 2025, 9:55pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329 "2025-10-16T21:55:52Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![ngu](https://avatars.discourse-cdn.com/v4/letter/n/58f4c7/32.png) [@ngu](https://community.geodynamics.org/u/ngu)\
**Post date:** [October 16, 2025, 9:55pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/1 "2025-10-16T21:55:53Z")

</div>

I am trying the reactive fluid tranport ans i have the following errors:  
void aspect::MaterialModel::ReactiveFluidTransport::melt\_fractions(const aspect::MaterialModel::MaterialModelInputs&, std::vector&, const aspect::MaterialModel::MaterialModelOutputs\*) const [with int dim = 3]  
The violated condition was:  
out != nullptr  
Additional information:  
The material model ‘ReactiveFluidTransport’ requires the material  
model outputs in order to compute melt fractions, but none were  
provided.  
I wish to find out if the i have made an error somewhere with this initial composition model:  
#include \<aspect/initial\_composition/lithosphere\_asthenosphere.h\>  
#include \<aspect/initial\_temperature/interface.h\>  
#include \<aspect/postprocess/interface.h\>  
#include \<aspect/utilities.h\>  
#include \<deal.II/base/signaling\_nan.h\>  
#include \<aspect/adiabatic\_conditions/interface.h\>  
#include \<aspect/material\_model/interface.h\>  
#include \<aspect/melt.h\>  
#include \<aspect/geometry\_model/interface.h\>

namespace aspect  
{  
namespace InitialComposition  
{

template   
LithosphereAsthenosphere::LithosphereAsthenosphere()  
:  
surface\_boundary\_id(5)  
{}

template   
void  
LithosphereAsthenosphere::initialize ()  
{

// Utilities::AsciiDataInitial::initialize(1);  
// Utilities::AsciiDataBoundary::initialize(1);  
surface\_boundary\_id = this-\>get\_geometry\_model().translate\_symbolic\_boundary\_name\_to\_id(“outer”);

```
std::set<types::boundary_id> surface_boundary_set;
surface_boundary_set.insert(surface_boundary_id);

  Utilities::AsciiDataBoundary<dim>::initialize(surface_boundary_set,1);

template <int dim>
double
LithosphereAsthenosphere<dim>::

```

initial\_composition (const Point &position,  
const unsigned int n\_comp) const  
{  
AssertThrow(this-\>introspection().compositional\_name\_exists(“porosity”),  
ExcMessage("The initial composition plugin lithosphere porosity’ did not find a "  
"compositional field called porosity’ to initialize. Please add a "  
“compositional field with this name.”));

```
const MaterialModel::MeltFractionModel<dim> *material_model =
    dynamic_cast<const MaterialModel::MeltFractionModel<dim>* > (&this->get_material_model());

AssertThrow(material_model != nullptr,
            ExcMessage("The used material model is not derived from the 'MeltFractionModel' class, "
                       "and therefore does not support computing equilibrium melt fractions. "
                       "This is incompatible with the porosity' "
                       "initial composition plugin, which needs to compute these melt fractions."));

const unsigned int n_fields = this->n_compositional_fields();
MaterialModel::MaterialModelInputs<dim> in(1, n_fields); // Use the class member to initialise size properly

in.position[0] = position;
in.temperature[0] = this->get_initial_temperature_manager().initial_temperature(position);
in.pressure[0] = this->get_adiabatic_conditions().pressure(position);
in.pressure_gradient[0] = 0.0;
in.velocity[0] = 0.0;
in.strain_rate[0] = SymmetricTensor<2,dim>();

const double depth = this->get_geometry_model().depth(position);
const unsigned int porosity_index = this->introspection().compositional_index_for_name("porosity");
const unsigned int peridotite_index = this->introspection().compositional_index_for_name("peridotite");
const unsigned int bound_fluid_index = this->introspection().compositional_index_for_name("bound_fluid");

double LAB_depth;

```

// LAB\_depth = ascii\_data\_lab.get\_data\_component(surface\_boundary\_id, position, 0);   
LAB\_depth = Utilities::AsciiDataBoundary::get\_data\_component(surface\_boundary\_id, position, 0);

```
  if (depth < moho && n_comp == 0)
    return 1.;
  else if (depth >= moho && depth < LAB_depth && n_comp == 1)
    return 1.;

  else if (depth > LAB_depth && n_comp == porosity_index)
    {
      in.composition[0][porosity_index] = 0.0;

      in.strain_rate[0] = SymmetricTensor<2,dim>();

      std::vector<double> melt_fraction(1);
      material_model->melt_fractions(in,melt_fraction);
      return melt_fraction[0];
    }

  else if (depth > LAB_depth && n_comp == peridotite_index)
    {
      in.composition[0][peridotite_index] = 1.0;

      in.strain_rate[0] = SymmetricTensor<2,dim>();

      std::vector<double> melt_fraction(1);
      material_model->melt_fractions(in,melt_fraction);
      return 1 - melt_fraction[0];
    }

  else
    return 0.;
}
template <int dim>
void
LithosphereAsthenosphere<dim>::declare_parameters(ParameterHandler &prm)
{
  prm.enter_subsection ("Initial composition model");
  {
    prm.enter_subsection("Lithosphere asthenosphere");
    {
      prm.declare_entry ("Moho", "30000.0",
                         Patterns::Double (0),
                         "Moho depth. Units: $m.");

         prm.enter_subsection("LAB file");
         {

```

// Utilities::AsciiDataBoundary::declare\_parameters(prm,  
Utilities::AsciiDataBase::declare\_parameters(prm,  
“$ASPECT\_SOURCE\_DIR/data/initial-temperature/ascii-data/”,  
“litho.africa.fishwick.txt”);  
// “litho.africa.200km.txt”);  
}  
prm.leave\_subsection();  
}  
prm.leave\_subsection();  
}  
prm.leave\_subsection();  
}

```
template <int dim>
void
LithosphereAsthenosphere<dim>::parse_parameters (ParameterHandler &prm)
{

```

// Utilities::AsciiDataBase::parse\_parameters(prm);

```
  prm.enter_subsection("Initial composition model");
  {
    prm.enter_subsection("Lithosphere asthenosphere");
    {
      moho = prm.get_double ("Moho");

```

// LAB\_isotherm = prm.get\_double (“LAB isotherm”);  
prm.enter\_subsection(“LAB file”);  
{  
Utilities::AsciiDataBase::parse\_parameters(prm);  
// ascii\_data\_lab.initialize\_simulator (this-\>get\_simulator());  
// ascii\_data\_lab.parse\_parameters(prm);  
}  
prm.leave\_subsection();

```
    }
    prm.leave_subsection();
  }
  prm.leave_subsection();
}

```

}  
}

// explicit instantiations  
namespace aspect  
{  
namespace InitialComposition  
{  
ASPECT\_REGISTER\_INITIAL\_COMPOSITION\_MODEL(LithosphereAsthenosphere,  
“lithosphere asthenosphere”,  
"Specify the composition in terms of an explicit formula. The format of these "  
"functions follows the syntax understood by the "  
“muparser library, see Section~\ref{sec:muparser-format}.”)  
}  
}

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [October 20, 2025, 8:26pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/2 "2025-10-20T20:26:03Z")

</div>

Hi @ngu,

Thank you for posting the question to the forum. I confess I’m having some trouble reading through the post as currently constructed.

Can you either reformat the post so that all of the code has consistent formatting (put it all under a pre-formatted text block or quotes) or post the .cc and .prm files as attachments?

Aside, I think the underlying error message is fairly clear - you are trying to get the initial melt fraction from the material model within the initial composition model, but you have not provided the material model outputs from the base model (required for the reactive fluid transport model).

Perhaps it would also be good to provide a high-level overview of what you are trying to achieve, as there may be a way to simplify the code a bit or use an alternative approach.

John

---

<div class="post-metadata">

**Author:** ![ngu](https://avatars.discourse-cdn.com/v4/letter/n/58f4c7/32.png) [@ngu](https://community.geodynamics.org/u/ngu)\
**Post date:** [October 23, 2025, 2:42pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/3 "2025-10-23T14:42:01Z")

</div>

Thank you for the reply. This is the error

void aspect::MaterialModel::ReactiveFluidTransport::melt\_fractions(const aspect::MaterialModel::MaterialModelInputs&, std::vector&, const aspect::MaterialModel::MaterialModelOutputs\*) const [with int dim = 3]  
The violated condition was:  
out != nullptr  
Additional information:  
The material model ‘ReactiveFluidTransport’ requires the material  
model outputs in order to compute melt fractions, but none were  
provided.

I wish to find out if there is a problem with this initial composition. I am using the reactive fluid transport material model.

```auto
#include <aspect/initial_composition/lithosphere_asthenosphere.h>
#include <aspect/initial_temperature/interface.h>
#include <aspect/postprocess/interface.h>
#include <aspect/utilities.h>
#include <deal.II/base/signaling_nan.h>
#include <aspect/adiabatic_conditions/interface.h>
#include <aspect/material_model/interface.h>
#include <aspect/melt.h>
#include <aspect/geometry_model/interface.h>

namespace aspect
{
  namespace InitialComposition
  {

   template <int dim>
   LithosphereAsthenosphere<dim>::LithosphereAsthenosphere()
   :
   surface_boundary_id(5)
   {}

   template <int dim>
   void
   LithosphereAsthenosphere<dim>::initialize ()
   {

// Utilities::AsciiDataInitial<dim>::initialize(1);
// Utilities::AsciiDataBoundary<dim>::initialize(1);
    surface_boundary_id = this->get_geometry_model().translate_symbolic_boundary_name_to_id("outer");

    std::set<types::boundary_id> surface_boundary_set;
    surface_boundary_set.insert(surface_boundary_id);

      Utilities::AsciiDataBoundary<dim>::initialize(surface_boundary_set,1);

    template <int dim>
    double
    LithosphereAsthenosphere<dim>::
initial_composition (const Point<dim> &position,
                     const unsigned int n_comp) const
{
    AssertThrow(this->introspection().compositional_name_exists("porosity"),
                ExcMessage("The initial composition plugin lithosphere porosity' did not find a "
                           "compositional field called porosity' to initialize. Please add a "
                           "compositional field with this name."));

    const MaterialModel::MeltFractionModel<dim> *material_model =
        dynamic_cast<const MaterialModel::MeltFractionModel<dim>* > (&this->get_material_model());

    AssertThrow(material_model != nullptr,
                ExcMessage("The used material model is not derived from the 'MeltFractionModel' class, "
                           "and therefore does not support computing equilibrium melt fractions. "
                           "This is incompatible with the porosity' "
                           "initial composition plugin, which needs to compute these melt fractions."));

    const unsigned int n_fields = this->n_compositional_fields();
    MaterialModel::MaterialModelInputs<dim> in(1, n_fields); // Use the class member to initialise size properly

    in.position[0] = position;
    in.temperature[0] = this->get_initial_temperature_manager().initial_temperature(position);
    in.pressure[0] = this->get_adiabatic_conditions().pressure(position);
    in.pressure_gradient[0] = 0.0;
    in.velocity[0] = 0.0;
    in.strain_rate[0] = SymmetricTensor<2,dim>();

    const double depth = this->get_geometry_model().depth(position);
    const unsigned int porosity_index = this->introspection().compositional_index_for_name("porosity");
    const unsigned int peridotite_index = this->introspection().compositional_index_for_name("peridotite");
    const unsigned int bound_fluid_index = this->introspection().compositional_index_for_name("bound_fluid");

    double LAB_depth;

// LAB_depth = ascii_data_lab.get_data_component(surface_boundary_id, position, 0);	
       	   LAB_depth = Utilities::AsciiDataBoundary<dim>::get_data_component(surface_boundary_id, position, 0);

      if (depth < moho && n_comp == 0)
        return 1.;
      else if (depth >= moho && depth < LAB_depth && n_comp == 1)
        return 1.;
    
      else if (depth > LAB_depth && n_comp == porosity_index)
        {
          in.composition[0][porosity_index] = 0.0;
    
          in.strain_rate[0] = SymmetricTensor<2,dim>();

          std::vector<double> melt_fraction(1);
          material_model->melt_fractions(in,melt_fraction);
          return melt_fraction[0];
        }

      else if (depth > LAB_depth && n_comp == peridotite_index)
        {
          in.composition[0][peridotite_index] = 1.0;

          in.strain_rate[0] = SymmetricTensor<2,dim>();
    
          std::vector<double> melt_fraction(1);
          material_model->melt_fractions(in,melt_fraction);
          return 1 - melt_fraction[0];
        }

      else
        return 0.;
    }
    template <int dim>
    void
    LithosphereAsthenosphere<dim>::declare_parameters(ParameterHandler &prm)
    {
      prm.enter_subsection ("Initial composition model");
      {
        prm.enter_subsection("Lithosphere asthenosphere");
        {
          prm.declare_entry ("Moho", "30000.0",
                             Patterns::Double (0),
                             "Moho depth. Units: $m.");

             prm.enter_subsection("LAB file");
             {
// Utilities::AsciiDataBoundary<dim>::declare_parameters(prm,
                 Utilities::AsciiDataBase<dim>::declare_parameters(prm,
                                                            "$ASPECT_SOURCE_DIR/data/initial-temperature/ascii-data/",
                                                             "litho.africa.fishwick.txt");
// "litho.africa.200km.txt");
             }
            prm.leave_subsection();
       	}
       	prm.leave_subsection();
	}
	prm.leave_subsection();
	}

    template <int dim>
    void
    LithosphereAsthenosphere<dim>::parse_parameters (ParameterHandler &prm)
    {
// Utilities::AsciiDataBase<dim>::parse_parameters(prm);

      prm.enter_subsection("Initial composition model");
      {
        prm.enter_subsection("Lithosphere asthenosphere");
        {
          moho = prm.get_double ("Moho");
// LAB_isotherm = prm.get_double ("LAB isotherm");
         prm.enter_subsection("LAB file");
         {
        Utilities::AsciiDataBase<dim>::parse_parameters(prm);
// ascii_data_lab.initialize_simulator (this->get_simulator());
// ascii_data_lab.parse_parameters(prm);
         }
         prm.leave_subsection();

        }
        prm.leave_subsection();
      }
      prm.leave_subsection();
    }

  }
}

// explicit instantiations
namespace aspect
{
  namespace InitialComposition
  {
    ASPECT_REGISTER_INITIAL_COMPOSITION_MODEL(LithosphereAsthenosphere,
                                              "lithosphere asthenosphere",
                                              "Specify the composition in terms of an explicit formula. The format of these "
                                              "functions follows the syntax understood by the "
                                              "muparser library, see Section~\\ref{sec:muparser-format}.")
  }
}

```

---

<div class="post-metadata">

**Author:** ![ngu](https://avatars.discourse-cdn.com/v4/letter/n/58f4c7/32.png) [@ngu](https://community.geodynamics.org/u/ngu)\
**Post date:** [October 23, 2025, 2:43pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/4 "2025-10-23T14:43:50Z")

</div>

Thanks for the reply John. I have provided a better version of the problem above

---

<div class="post-metadata">

**Author:** ![daniel.douglas](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/daniel.douglas/32/2126_2.png) [@daniel.douglas](https://community.geodynamics.org/u/daniel.douglas)\
**Post date:** [October 27, 2025, 5:05pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/5 "2025-10-27T17:05:31Z")

</div>

Hi @ngu

Thanks for posting the cleaned up version of your plugin! I agree with John that the error message that is being given is informative. There are several places in your plugin where you are calling the `melt_fractions` function without providing the necessary arguments. You can see exactly what arguments `ReactiveFluidTransport.melt_fractions` requires [here](https://github.com/geodynamics/aspect/blob/main/source/material_model/reactive_fluid_transport.cc#L82), and see that it requires `MaterialModelInputs`, as well as `MaterialModelOutputs`.

Cheers,

Daniel

---

<div class="post-metadata">

**Author:** ![ngu](https://avatars.discourse-cdn.com/v4/letter/n/58f4c7/32.png) [@ngu](https://community.geodynamics.org/u/ngu)\
**Post date:** [October 27, 2025, 5:42pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/6 "2025-10-27T17:42:16Z")

</div>

Thanks again John. I am trying to generate melt within the asthernsophere using the reactive fluid transport material model with 5 compositions, crust, mantle lithosphere, porosity, peridotite and bound fluid

---

<div class="post-metadata">

**Author:** ![ngu](https://avatars.discourse-cdn.com/v4/letter/n/58f4c7/32.png) [@ngu](https://community.geodynamics.org/u/ngu)\
**Post date:** [October 27, 2025, 5:43pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/7 "2025-10-27T17:43:12Z")

</div>

Thanks for the hint Daniel. Will look at it again and see if i can fix it

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [October 27, 2025, 7:54pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/8 "2025-10-27T19:54:10Z")

</div>

@daniel.douglas - Thanks for weighing in and outlining what inputs are required for that function call.

@ngu - I think this plugin could be simplified quite a bit by not specifying an initial porosity, and rather having that value dynamically calculated during the first time step. In that scenario, you would not need to interface with the material model. Is there are particular reason why you would like to have the porosity (melt fraction in this case) determined in the composition initialization stage?

---

<div class="post-metadata">

**Author:** ![ngu](https://avatars.discourse-cdn.com/v4/letter/n/58f4c7/32.png) [@ngu](https://community.geodynamics.org/u/ngu)\
**Post date:** [November 7, 2025, 9:14pm UTC](https://community.geodynamics.org/t/problems-with-initial-composition-model-using-reactive-fluid-transport/4329/9 "2025-11-07T21:14:07Z")

</div>

Thanks very much John  
I tried to do that without using the initial porosity but I got “illegal instruction (core dumped)” as the error. Maybe am doing something wrong. Can you please help me direct how i can do that without the initial porosity?