Create reactive fluid transport material model #5245

naliboff · 2023-07-10T03:37:03Z

A draft material model, which composites on top of existing material models for solid deformation to also simulate reactive fluid transport. This is based on the contribution from #5214 and discussions with @jdannberg, @danieldouglas92, and @ryanstoner1.

The material is a work in progress and the PR was created to facilitate discussion. At present, the only option is to simulate non-reactive (insoluble) fluid flow, but the idea is to maintain a flexible structure and allow users to use add additional petrologic models of varying complexity.

I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

I have tested my new feature locally to ensure it is correct.
I have created a testcase for the new feature/benchmark in the tests/ directory.
I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

jdannberg

Thank you! I have looked through your PR and had some smaller comments, and I also think it would be good if you could reorder the functions. Usually, all material models have the same order of the functions (for example, the last 4 are always evaluate(), declare_parameters(), parse_parameters(), and create_additional_named_outputs()). Have a look at the other melt material models for the right order.

Maybe it would also be good if someone else could have a look (just to comment on if the general idea of this makes sense since it duplicates some code).

include/aspect/material_model/reactive_fluid_transport.h

source/material_model/reactive_fluid_transport.cc

jdannberg · 2023-07-10T05:30:08Z

source/material_model/reactive_fluid_transport.cc

+          // the melt (fluid) fraction at any point is thus equal to the
+          // porosity, which is determined by the assigned initial porosity,
+          // fluid boundary conditions, and fluid transport through the model.


I am not sure that this is true. If there is bound water, the total fluid fraction would not equal the porosity.
The free fluid fraction is equal to the porosity.

Good point, I've modified the documentation and code accordingly.

source/material_model/reactive_fluid_transport.cc

jdannberg · 2023-07-10T05:33:06Z

tests/reactive_fluid_transport.prm

@@ -0,0 +1,155 @@
+# This is a 1D test case for an insoluble fluid moving upwards in a
+# narrow channel through porous coupled two-phase flow.


The tester says: "test reactive_fluid_transport.prm is missing a folder with reference data"

tests/reactive_fluid_transport.prm

jdannberg · 2023-07-10T05:34:48Z

tests/reactive_fluid_transport.prm

+# This is a 1D test case for an insoluble fluid moving upwards in a
+# narrow channel through porous coupled two-phase flow.


I think this would be more correct

Suggested change

# This is a 1D test case for an insoluble fluid moving upwards in a

# narrow channel through porous coupled two-phase flow.

# This is a 1D test case for a fluid moving upwards in a narrow channel

#through porous coupled two-phase flow without reacting with the solid rock.

jdannberg · 2023-07-10T05:35:34Z

tests/reactive_fluid_transport.prm

+# We use a custom material model that implements the layers with
+# different solubility.


This is outdated

jdannberg · 2023-07-10T05:35:50Z

tests/reactive_fluid_transport.prm

+    set Reference fluid viscosity = 10
+    set Reference permeability = 2.5e-6
+    set Exponential fluid weakening factor = 0
+    set Fluid reaction time scale for operator splitting = 5e4


But operator splitting is not used?

Good catch, removed

naliboff · 2023-07-10T15:52:44Z

@jdannberg - Thanks for the initial comments and suggestions, which I've addressed aside from the reordering the functions. I'll do the reordering next, finish the test case, and then add in a structure for specifying what reaction model to use. After that, a look from a second reviewer would be ideal.

naliboff · 2023-07-10T19:09:31Z

@jdannberg @ryanstoner1 @danieldouglas92 - The functions have been re-ordered, a switch case is available for implementing different reaction schemes (current option is no reaction) and the test results have been added. It is probably ready for another round of reviews?

@bobmyhill - Would you be willing to take a look?

naliboff · 2023-07-10T19:10:12Z

/rebuild

naliboff · 2023-07-10T19:13:25Z

Topics for discussion when reviewing:

What features should be present before this PR is merged (level of complexity)
Will the current structure facilitate future expansion?
Should the test case be modified to be based on a benchmark for simple Darcy flow?

ryanstoner1

@naliboff Some comments about settling on the semantics. Take them of leave them. Perhaps we can discuss some when we talk. Also, it might be good to go over how to best make this extensible if new users have their own reaction formulation.

ryanstoner1 · 2023-07-11T01:15:58Z

source/material_model/reactive_fluid_transport.cc

+              case no_reaction:
+              {
+                // A fluid-rock reaction model where no reactions occur.
+                // The melt (fluid) fraction at any point is equal


Perhaps specify this to make a bit more transparent?

Suggested change

// The melt (fluid) fraction at any point is equal

// The melt (fluid) fraction (volume percent) at any point is equal

ryanstoner1 · 2023-07-11T01:19:03Z

source/material_model/reactive_fluid_transport.cc

+      // Modify the viscosity from the base model based on the presence of fluid.
+      const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
+


Suggested change

// Modify the viscosity from the base model based on the presence of fluid.

const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");

const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");

// Modify the viscosity from the base model based on the presence of fluid.

ryanstoner1 · 2023-07-11T01:21:02Z

source/material_model/reactive_fluid_transport.cc

+                  const double phi_0 = 0.05;
+                  porosity = std::max(porosity,1e-8);


Comment could be nice here.

ryanstoner1 · 2023-07-11T01:26:50Z

source/material_model/reactive_fluid_transport.cc

+        {
+          prm.declare_entry("Base model","visco plastic",
+                            Patterns::Selection(MaterialModel::get_valid_model_names_pattern<dim>()),
+                            "The name of a material model that will be modified by the "


Suggested change

"The name of a material model that will be modified by the "

"The name of a material model incorporating the "

Clarity.

gassmoeller

I talked to @jdannberg and the general approach makes sense. Go ahead and address the comments.

jdannberg

Thank you! The structure looks good now, but I think the implementation requires some changes. Most importantly, think about if you actually want no reactions, a zero solubility, or both.

jdannberg · 2023-07-15T03:29:27Z

doc/modules/changes/20230714_naliboff_b

@@ -0,0 +1,8 @@
+New: There is now a new material ('reactive fluid 
+transport') that is designed to advect fluids and 
+compute luid release and absorption based on 


Suggested change

compute luid release and absorption based on

compute fluid release and absorption based on

jdannberg · 2023-07-15T03:29:37Z

doc/modules/changes/20230714_naliboff_b

@@ -0,0 +1,8 @@
+New: There is now a new material ('reactive fluid 


Suggested change

New: There is now a new material ('reactive fluid

New: There is now a new material model ('reactive fluid

jdannberg · 2023-07-15T03:41:09Z

doc/modules/changes/20230714_naliboff_b

+transport') that is designed to advect fluids and 
+compute luid release and absorption based on 
+different models for fluid-rock interaction. The
+that properties of the solid are taken from a 


Suggested change

that properties of the solid are taken from a

properties of the solid are taken from a

jdannberg · 2023-07-15T03:44:10Z

include/aspect/material_model/reactive_fluid_transport.h

+        double fluid_reaction_time_scale;
+
+        /**
+         * Enumeration for selecting which type of scheme to use for reactions between.


Suggested change

* Enumeration for selecting which type of scheme to use for reactions between.

* Enumeration for selecting which type of scheme to use for reactions between

jdannberg · 2023-07-15T03:45:28Z

source/material_model/reactive_fluid_transport.cc

+              {
+                // A fluid-rock reaction model where no reactions occur.
+                // The melt (fluid) fraction (volume percent) at any point is
+                // equal to the sume of the bound fluid content and porosity,


Suggested change

// equal to the sume of the bound fluid content and porosity,

// equal to the sum of the bound fluid content and porosity,

jdannberg · 2023-07-15T04:17:24Z

include/aspect/material_model/reactive_fluid_transport.h

+         */
+        enum ReactionScheme
+        {
+          no_reaction


jdannberg · 2023-07-15T04:22:03Z

source/material_model/reactive_fluid_transport.cc

+                             "computed. If the model does not use operator splitting, this parameter is not used. "
+                             "Units: yr or s, depending on the ``Use years "
+                             "in output instead of seconds'' parameter.");
+          prm.declare_entry ("Fluid solid reaction scheme", "no reaction",


I think I would prefer this with a dash

Suggested change

prm.declare_entry ("Fluid solid reaction scheme", "no reaction",

prm.declare_entry ("Fluid-solid reaction scheme", "no reaction",

jdannberg · 2023-07-15T04:24:51Z

tests/reactive_fluid_transport_no_reaction/screen-output

+
+
+   Postprocessing:
+     Compositions min/max/mass:   -5.788e-07/0.11/5159 // 0/0/0


Are you sure this is working as intended? Shouldn't the fluid be released, i.e. the porosity should increase?

Yep, good catch! The initial assignment of values (porosity, bound_fluid) was reversed.

jdannberg · 2023-07-15T04:29:03Z

source/material_model/reactive_fluid_transport.cc

+          for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
+            for (unsigned int c=0; c<in.composition[q].size(); ++c)
+              {
+                double porosity_change = eq_free_fluid_fractions[q] - in.composition[q][porosity_idx];
+                // do not allow negative porosity
+                if (in.composition[q][porosity_idx] + porosity_change < 0)
+                  porosity_change = -in.composition[q][porosity_idx];
+
+                if (c == bound_fluid_idx && this->get_timestep_number() > 0)
+                  reaction_rate_out->reaction_rates[q][c] = - porosity_change / fluid_reaction_time_scale;
+                else if (c == porosity_idx && this->get_timestep_number() > 0)
+                  reaction_rate_out->reaction_rates[q][c] = porosity_change / fluid_reaction_time_scale;
+                else
+                  reaction_rate_out->reaction_rates[q][c] = 0.0;
+              }


This only works if you use operator splitting in your model. So (1), please document that, and (2) change your test to actually use operator splitting (because at the moment, it does not test this implementation because it does not use the reaction rates).

Can you elaborate on what part of the code needs additional documentation for the operator splitting? There is some documentation at the top on this above the following if statement (if (this->get_parameters().use_operator_splitting && reaction_rate_out != nullptr)), but happy to more in an additional location.

The test has been updated

Just having the assert is fine!

jdannberg · 2023-07-15T04:29:31Z

tests/reactive_fluid_transport_no_reaction.prm

+set End time                               = 1e1
+
+
+set Use operator splitting                     = false


You need to use operator splitting for your reactions to work.

naliboff · 2023-10-07T20:04:42Z

@jdannberg - Thanks for the thorough review! I've hopefully addressed most of the comments, and I think this is ready for another review.

jdannberg

I only have a few minor comments, otherwise this looks good!

jdannberg · 2023-10-07T20:52:21Z

include/aspect/material_model/reactive_fluid_transport.h

+         * model is thus assuming that the bound water fraction is a
+         * volume fraction (i.e., since porosity is always a volume
+         * fraction). This latter assumption also requires the selected
+         * base model is incompressible, as otherwise the advection


Suggested change

* base model is incompressible, as otherwise the advection

* base model to be incompressible, as otherwise the advection

jdannberg · 2023-10-07T20:54:22Z

tests/reactive_fluid_transport_no_reaction.prm

+# there is no need to use operator splitting.
+set Use operator splitting                     = false
+
+# Mellt transport is turned on to ensure it works


Suggested change

# Mellt transport is turned on to ensure it works

# Melt transport is turned on to ensure it works

jdannberg · 2023-10-07T21:00:45Z

tests/reactive_fluid_transport_no_reaction.prm

+    set Reference fluid viscosity = 10
+    set Reference permeability = 2.5e-6
+    set Exponential fluid weakening factor = 0
+    set Fluid-solid reaction scheme = zero solubility


Shouldn't this be no reaction?
You switch the operator splitting off, so technically, this should give you the correct result, but I find it to be confusing.

I see it actually does not work because
Exception 'ExcMessage("The Fluid-reaction scheme zero solubility must be used with operator splitting.")'
:-)

jdannberg · 2023-10-07T21:02:21Z

tests/reactive_fluid_transport_zero_solubility.prm

@@ -0,0 +1,161 @@
+


Can you remove this empty line?

jdannberg · 2023-10-07T21:05:59Z

tests/reactive_fluid_transport_no_reaction.prm

+# This model tests the reactive fluid transport material
+# 'no reaction' scenario to ensure that there is no reaction
+# between the fluid and solid phases. The values of
+# porosity (0) and the bound fluid (0.01) are initially
+# spatially uniform. The absence of reactions between
+# the two phases and near zero velocities (free-slip on
+# all sides) produces no change in their values through time.


This is very similar to the zero solubility test below. Can you include the other test here with
include $ASPECT_SOURCE_DIR/tests/reactive_fluid_transport_zero_solubility.prm
and the only list the parameters that are different? This way you do not have to duplicate everything.

jdannberg · 2023-10-07T21:09:05Z

tests/reactive_fluid_transport_zero_solubility/screen-output

+     Compositions min/max/mass:   0.01/0.01/4e+04 // 0/0/0
+     Melt fraction min/total/max: 0.01, 4e+04, 0.01


This shows nicely that it worked!

jdannberg · 2023-10-07T21:18:27Z

source/material_model/reactive_fluid_transport.cc

+          for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
+            for (unsigned int c=0; c<in.composition[q].size(); ++c)
+              {
+                double porosity_change = eq_free_fluid_fractions[q] - in.composition[q][porosity_idx];
+                // do not allow negative porosity
+                if (in.composition[q][porosity_idx] + porosity_change < 0)
+                  porosity_change = -in.composition[q][porosity_idx];
+
+                if (c == bound_fluid_idx && this->get_timestep_number() > 0)
+                  reaction_rate_out->reaction_rates[q][c] = - porosity_change / fluid_reaction_time_scale;
+                else if (c == porosity_idx && this->get_timestep_number() > 0)
+                  reaction_rate_out->reaction_rates[q][c] = porosity_change / fluid_reaction_time_scale;
+                else
+                  reaction_rate_out->reaction_rates[q][c] = 0.0;
+              }


Just having the assert is fine!

jdannberg · 2023-10-07T21:41:23Z

source/material_model/reactive_fluid_transport.cc

+            {
+              AssertThrow(fluid_reaction_time_scale >= this->get_parameters().reaction_time_step,
+                          ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
+                                     + " in the operator splitting scheme is too large to compute melting rates! "


Suggested change

+ " in the operator splitting scheme is too large to compute melting rates! "

+ " in the operator splitting scheme is too large to compute fluid release rates! "

jdannberg · 2023-10-07T21:41:42Z

source/material_model/reactive_fluid_transport.cc

+              AssertThrow(fluid_reaction_time_scale >= this->get_parameters().reaction_time_step,
+                          ExcMessage("The reaction time step " + Utilities::to_string(this->get_parameters().reaction_time_step)
+                                     + " in the operator splitting scheme is too large to compute melting rates! "
+                                     "You have to choose it in such a way that it is smaller than the 'Melting time scale for "


Suggested change

"You have to choose it in such a way that it is smaller than the 'Melting time scale for "

"You have to choose it in such a way that it is smaller than the 'Fluid reaction time scale for "

jdannberg · 2023-10-07T21:44:30Z

source/material_model/reactive_fluid_transport.cc

+              AssertThrow(this->get_material_model().is_compressible() == false,
+                          ExcMessage("The Fluid-reaction scheme zero solubility must be used with an incompressible base model."));


This does not work well for all models, since at this point we have not parsed the parameters from the base model. So if the base model has a compressibility that depends on an input parameter, at this point in the code we will not know about it. It would be better to move this after parsing the parameters from the base model or ito the initialize() function.

naliboff · 2023-10-07T23:03:25Z

@jdannberg - Thanks for the last round of edits! I addressed all of the comments, but there is still a bug I am tracking down so not ready to merge yet.

naliboff · 2023-12-08T19:04:15Z

@jdannberg - I think this is ready to merged!

jdannberg · 2024-01-30T18:38:35Z

This looks good to me now! The only other comment I have is that this also does not create the additional material model outputs the base model has (which I fixed for the prescribed viscosity and depth dependent model here #5547).

But that could be addressed in a follow-up PR.

And sorry this took me so long to look at!

jdannberg reviewed Jul 10, 2023

View reviewed changes

naliboff force-pushed the create_reactive_fluid_transport_material_model branch from 3d5b8ae to 1fdafa8 Compare July 10, 2023 15:49

naliboff added the ready to test label Jul 10, 2023

ryanstoner1 reviewed Jul 11, 2023

View reviewed changes

gassmoeller reviewed Jul 13, 2023

View reviewed changes

naliboff changed the title ~~[WIP] Create reactive fluid transport material model~~ Create reactive fluid transport material model Jul 14, 2023

naliboff added 7 commits July 14, 2023 11:59

add reactive fluid transport material model

7ad2517

add test case

2383272

address comments

2270160

reorder functions, add parameter case switch

2c13c93

update test file and results

acdd9f1

fix parameter documentation

b8843c5

address comments

1318d47

naliboff force-pushed the create_reactive_fluid_transport_material_model branch from 0147f46 to 4a69158 Compare July 14, 2023 18:10

naliboff added 3 commits July 14, 2023 14:10

fix compiler error

4dd7f43

add change log

f688eb0

fix another compiler error

5e5bd1c

naliboff force-pushed the create_reactive_fluid_transport_material_model branch from 4a69158 to 5e5bd1c Compare July 14, 2023 20:10

naliboff mentioned this pull request Jul 14, 2023

Add Tian Parametrization to Reactive Fluid Transport #5313

Merged

5 tasks

jdannberg requested changes Jul 15, 2023

View reviewed changes

rework reactive transport zero solubility option

c4c283e

naliboff force-pushed the create_reactive_fluid_transport_material_model branch from bc91b01 to 2d54a9f Compare October 6, 2023 23:23

add reactive transport no reaction option

ce8be6f

naliboff force-pushed the create_reactive_fluid_transport_material_model branch from 2d54a9f to ce8be6f Compare October 6, 2023 23:40

update test

71a5a05

naliboff force-pushed the create_reactive_fluid_transport_material_model branch from 60310f6 to 71a5a05 Compare October 6, 2023 23:58

naliboff added 2 commits October 7, 2023 11:59

add asserts

7c2abdb

cleanup reactive transport material

1cb40d3

jdannberg reviewed Oct 7, 2023

View reviewed changes

address comments

b90fb04

naliboff added 2 commits October 7, 2023 17:17

fix a conditional statement

2c03879

update tests

ae7eb29

jdannberg approved these changes Jan 30, 2024

View reviewed changes

jdannberg merged commit 3f38e93 into geodynamics:main Jan 31, 2024
8 checks passed

		@@ -0,0 +1,155 @@
		# This is a 1D test case for an insoluble fluid moving upwards in a
		# narrow channel through porous coupled two-phase flow.

		# We use a custom material model that implements the layers with
		# different solubility.

	// The melt (fluid) fraction at any point is equal
	// The melt (fluid) fraction (volume percent) at any point is equal

		// Modify the viscosity from the base model based on the presence of fluid.
		const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");

		const double phi_0 = 0.05;
		porosity = std::max(porosity,1e-8);

	"The name of a material model that will be modified by the "
	"The name of a material model incorporating the "

	compute luid release and absorption based on
	compute fluid release and absorption based on

		@@ -0,0 +1,8 @@
		New: There is now a new material ('reactive fluid

	New: There is now a new material ('reactive fluid
	New: There is now a new material model ('reactive fluid

	that properties of the solid are taken from a
	properties of the solid are taken from a

	* Enumeration for selecting which type of scheme to use for reactions between.
	* Enumeration for selecting which type of scheme to use for reactions between

	// equal to the sume of the bound fluid content and porosity,
	// equal to the sum of the bound fluid content and porosity,

	prm.declare_entry ("Fluid solid reaction scheme", "no reaction",
	prm.declare_entry ("Fluid-solid reaction scheme", "no reaction",



		Postprocessing:
		Compositions min/max/mass: -5.788e-07/0.11/5159 // 0/0/0

	* base model is incompressible, as otherwise the advection
	* base model to be incompressible, as otherwise the advection

	# Mellt transport is turned on to ensure it works
	# Melt transport is turned on to ensure it works

		Compositions min/max/mass: 0.01/0.01/4e+04 // 0/0/0
		Melt fraction min/total/max: 0.01, 4e+04, 0.01

	+ " in the operator splitting scheme is too large to compute melting rates! "
	+ " in the operator splitting scheme is too large to compute fluid release rates! "

	"You have to choose it in such a way that it is smaller than the 'Melting time scale for "
	"You have to choose it in such a way that it is smaller than the 'Fluid reaction time scale for "

		AssertThrow(this->get_material_model().is_compressible() == false,
		ExcMessage("The Fluid-reaction scheme zero solubility must be used with an incompressible base model."));

Create reactive fluid transport material model #5245

Create reactive fluid transport material model #5245

Conversation

naliboff commented Jul 10, 2023 • edited

For new features/models or changes of existing features:

jdannberg left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff commented Jul 10, 2023

naliboff commented Jul 10, 2023

naliboff commented Jul 10, 2023

naliboff commented Jul 10, 2023

ryanstoner1 left a comment

Choose a reason for hiding this comment

ryanstoner1 Jul 11, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gassmoeller left a comment

Choose a reason for hiding this comment

jdannberg left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff commented Oct 7, 2023

jdannberg left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff commented Oct 7, 2023

naliboff commented Dec 8, 2023

jdannberg commented Jan 30, 2024

naliboff commented Jul 10, 2023 •

edited

ryanstoner1 Jul 11, 2023 •

edited