@@ -86,25 +86,44 @@ namespace aspect
8686 std::vector<double > &melt_fractions,
8787 const MaterialModel::MaterialModelOutputs<dim> *out) const
8888 {
89- const std::shared_ptr< const MeltOutputs<dim>> fluid_out = out-> template get_additional_output_object <MeltOutputs<dim>>();
90- if (fluid_out != nullptr && in. requests_property (MaterialProperties::additional_outputs) )
89+ // the Katz 2003 model does not need information about fluid content
90+ if (fluid_solid_reaction_scheme == katz2003 )
9191 {
92+ for (unsigned int q=0 ; q<in.n_evaluation_points (); ++q)
93+ melt_fractions[q] = katz2003_model.melt_fraction (in.temperature [q],
94+ this ->get_adiabatic_conditions ().pressure (in.position [q]));
95+ }
96+ else
97+ {
98+ Assert (out != nullptr ,
99+ ExcMessage (" The material model 'ReactiveFluidTransport' requires the material model "
100+ " outputs in order to compute melt fractions, but none were provided." ));
101+ Assert (out->template has_additional_output_object <MeltOutputs<dim>>(),
102+ ExcMessage (" The material model 'ReactiveFluidTransport' requires melt material "
103+ " outputs to compute the melt fractions, but none were provided." ));
104+
105+ const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object <MeltOutputs<dim>>();
106+
107+ const unsigned int porosity_idx = this ->introspection ().compositional_index_for_name (" porosity" );
108+ const unsigned int bound_fluid_idx = (fluid_solid_reaction_scheme == zero_solubility)
109+ ?
110+ this ->introspection ().compositional_index_for_name (" bound_fluid" )
111+ :
112+ numbers::invalid_unsigned_int;
113+
92114 for (unsigned int q=0 ; q<in.n_evaluation_points (); ++q)
93115 {
94- const unsigned int porosity_idx = this ->introspection ().compositional_index_for_name (" porosity" );
116+ // This function outputs the mass fraction of equilibrium free water at each quadrature point. However,
117+ // the free water is a volume fraction, while the bound water is a mass fraction. Convert the free water
118+ // to a mass fraction so that we correctly compare the two values.
119+ const double volume_fraction_porosity = in.composition [q][porosity_idx];
120+ const double bulk_density = compute_bulk_density (volume_fraction_porosity, out->densities [q], fluid_out->fluid_densities [q]);
121+ const double mass_fraction_porosity = compute_mass_fraction (volume_fraction_porosity, fluid_out->fluid_densities [q], bulk_density);
122+
95123 switch (fluid_solid_reaction_scheme)
96124 {
97125 case no_reaction:
98126 {
99- Assert (out != nullptr ,
100- ExcMessage (" The MaterialModelOutputs provided to the No Reactions reaction model is a nullptr." ));
101- const double volume_fraction_porosity = in.composition [q][porosity_idx];
102- // Melt fractions outputs the mass fraction of equilibrium free water at each quadrature point. However,
103- // the free water is a volume fraction, while the bound water is a mass fraction. Convert the free water
104- // to a mass fraction so that we correctly compare the two values.
105- const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object <MeltOutputs<dim>>();
106- const double bulk_density = compute_bulk_density (volume_fraction_porosity, out->densities [q], fluid_out->fluid_densities [q]);
107- const double mass_fraction_porosity = compute_mass_fraction (volume_fraction_porosity, fluid_out->fluid_densities [q], bulk_density);
108127 // No reactions occur between the solid and fluid phases,
109128 // and the fluid mass fraction (stored in the melt_fractions
110129 // vector) is equal to the mass fraction of the porosity.
@@ -116,39 +135,14 @@ namespace aspect
116135 // The fluid volume fraction in equilibrium with the solid
117136 // at any point (stored in the melt_fractions vector) is
118137 // equal to the sum of the bound fluid content and porosity.
119- Assert (out != nullptr ,
120- ExcMessage (" The MaterialModelOutputs provided to the Zero Solubility reaction model is a nullptr." ));
121- const double volume_fraction_porosity = in.composition [q][porosity_idx];
122- const unsigned int bound_fluid_idx = this ->introspection ().compositional_index_for_name (" bound_fluid" );
123- // Melt fractions outputs the mass fraction of equilibrium free water at each quadrature point. However,
124- // the free water is a volume fraction, while the bound water is a mass fraction. Convert the free water
125- // to a mass fraction so that we correctly compare the two values.
126- const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object <MeltOutputs<dim>>();
127- const double bulk_density = compute_bulk_density (volume_fraction_porosity, out->densities [q], fluid_out->fluid_densities [q]);
128- const double mass_fraction_porosity = compute_mass_fraction (volume_fraction_porosity, fluid_out->fluid_densities [q], bulk_density);
129138 melt_fractions[q] = in.composition [q][bound_fluid_idx] + mass_fraction_porosity;
130139 break ;
131140 }
132141 case tian_approximation:
133142 {
134- Assert (out != nullptr ,
135- ExcMessage (" The MaterialModelOutputs provided to the Tian 2019 reaction model is a nullptr." ));
136- // The tian2019_model melt_fraction function requires the mass fraction of free water at each quadrature
137- // point as an input. However, the free water is stored in the porosity field as a volume fraction.
138- // We therefore convert the free water to a mass fraction so that we have the correct input.
139- const double volume_fraction_porosity = in.composition [q][porosity_idx];
140- const std::shared_ptr<const MeltOutputs<dim>> fluid_out = out->template get_additional_output_object <MeltOutputs<dim>>();
141- const double bulk_density = compute_bulk_density (volume_fraction_porosity, out->densities [q], fluid_out->fluid_densities [q]);;
142- const double mass_fraction_porosity = compute_mass_fraction (volume_fraction_porosity, fluid_out->fluid_densities [q], bulk_density);
143143 melt_fractions[q] = tian2019_model.melt_fraction (in, mass_fraction_porosity, q);
144144 break ;
145145 }
146- case katz2003:
147- {
148- melt_fractions[q] = katz2003_model.melt_fraction (in.temperature [q],
149- this ->get_adiabatic_conditions ().pressure (in.position [q]));
150- break ;
151- }
152146 default :
153147 {
154148 AssertThrow (false , ExcNotImplemented ());
@@ -236,6 +230,10 @@ namespace aspect
236230 if (this ->get_parameters ().use_operator_splitting && reaction_rate_out != nullptr
237231 && in.requests_property (MaterialProperties::reaction_rates))
238232 {
233+ Assert (fluid_out != nullptr ,
234+ ExcMessage (" The material model 'ReactiveFluidTransport' requires melt material "
235+ " outputs to compute reaction rates, but none were provided." ));
236+
239237 // Fill reaction rate outputs if the model uses operator splitting.
240238 // Specifically, change the porosity (representing the amount of free fluid)
241239 // based on the water solubility and the fluid content.
0 commit comments