Added reaction and melt output calculations to reactive fluid flow model #5795

Grant-Block · 2024-06-05T19:34:05Z

Pull Request Checklist. Please read and check each box with an X. Delete any part not applicable. Ask on the forum if you need help with any step.

This uses the KatzMantleMelting module to calculate reaction and melt terms as an option in the reactive_fluid_transport model. A test is included copying a MeltSimple test, but using a visco plastic solid material model.

Before your first pull request:

[ X] I have read the guidelines in our CONTRIBUTING.md document.

For all pull requests:

[ X] I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

[ X] I have tested my new feature locally to ensure it is correct.
[ X] 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. TODO

Grant-Block · 2024-06-05T19:35:17Z

Hi @jdannberg @naliboff, this PR builds off my other one, so it probably would be best to make sure the other one is wrapped up before doing too much more on this one. I just wanted to get this up!

jdannberg · 2024-06-05T21:40:38Z

Yes, you'll need to rebase to main before we can look at the changes (otherwise it's hard to tell which changes come from this PR and which ones are from the other one).

jdannberg

I had some comments on the general structure, and I'll look at the details once we have merged the other PR and this one is rebased.

jdannberg · 2024-06-05T22:22:41Z

source/material_model/reactive_fluid_transport.cc

@@ -377,6 +392,10 @@ namespace aspect
                             "use polynomials to describe hydration and dehydration reactions for four different "
                             "rock compositions as defined in Tian et al., 2019, or the Katz et. al. 2003 mantle "
                             "melting model.");


This material model now has a lot of parameters that are not being used when the katz2003 model is selected. The declaration will be overwritten by the ones in ReactionModel::Katz2003MantleMelting<dim>::declare_parameters(prm). For example, the default value of the "Reference permeability" will not be used and instead be overwritten by the one from the reaction model. There are two potential solutions

(1) If the computation of the melt outputs is similar, you remove the original parameters here, always call katz2003_model.calculate_fluid_outputs(in, out, reference_T), and make any small modifications necessary afterwards. That would be nice because it would also reduce the lines of code.
(2) If that does not work well, you could move the declare and parse parameters function of the ReactionModel::Katz2003MantleMelting into its own section in the parameter file, something like

prm.enter_subsection("Reactive Fluid Transport Model"); { prm.enter_subsection("Katz 2003 model"); {

This way, they would be separate.

I think I'm leaning towards the second solution because the Katz2003 model might not be the only melting model ReactiveFluidTransport has, I think I'll end up adding other melting models after the hackathon. Then it might make more sense for each melt model to have its own parameter subsection. Does that sound ok?

jdannberg · 2024-06-05T22:52:32Z

tests/viscoplastic_melt_blob_freezing/solution/solution-00002.pvtu

When you add tests with graphical output, change the output format to gnuplot rather than vtu and then use the gnuplot file. The reason is that vtu files are binary, so any small change in any number will change the whole file. With the gnuplot format, we can compare the inidividual numbers with numdiff.

…odel

Grant-Block · 2024-06-06T20:03:04Z

Hi @jdannberg and @naliboff, I rebased to main and squashed my commits here, so this PR is ready for comments. Thank you!

jdannberg

I only have a few small changes to suggest, otherwise this looks good to me! And have you checked that the output of the test is the same as for the rising_melt_blob_freezing test? Because that's what we would expect, right?

jdannberg · 2024-06-06T23:44:24Z

source/material_model/reactive_fluid_transport.cc

+                             "The reference temperature $T_0$. The reference temperature is used "
+                             "in both the density and viscosity formulas. Units: \\si{\\kelvin}.");


Suggested change

"The reference temperature $T_0$. The reference temperature is used "

"in both the density and viscosity formulas. Units: \\si{\\kelvin}.");

"The reference temperature $T_0$ for the katz2003 reaction model. The reference temperature is used "

"in both the density and viscosity formulas of this model. Units: \\si{\\kelvin}.");

jdannberg · 2024-06-06T23:50:48Z

tests/viscoplastic_melt_blob_freezing.prm

This is a copy of rising_melt_blob_freezing.prm, right? Instead of copying the test, can you include it by adding the line

include $ASPECT_SOURCE_DIR/tests/rising_melt_blob_freezing.prm

and then only including the lines that you changed compared to that test (this way, if something changes, we only have to update it in one test).

And can you also update the documentation? (Something like: "This is like the rising_melt_blob_freezing test, except that it uses the reactive fluid transport material model.")

Oh, and you still need to check why the tests are failing, and if those changes are only related to what you changed in the input file (the changes in shear_heating_with_melt looked tiny, but I couldn't tell how much has changed in viscoplastic_melt_fraction, if that was just related to adding the other field) and update the test results.

Grant-Block · 2024-06-07T00:01:59Z

Hi @jdannberg, thanks for your comments! I'll push the changes up to the PR after doing any changes @naliboff might have, so I don't make too many more commits! My new test has identical results to rising_melt_blob_freezing, except for the solid viscosity, which is what I would expect from using a different solid material in this case.

naliboff

@Grant-Block - Nice work!

Only two minor suggestions and a few self-reminder comments for future refactoring.

If possible, can you alter the test in this PR to give the same solid viscosity as the existing rising_melt_blob_freezing test (this way the tests are identical)? If not, can you make a note of the difference and why it is expected in the test PRM file?

Otherwise, this looks basically ready to go.

naliboff · 2024-06-06T21:16:01Z

source/material_model/reactive_fluid_transport.cc

        {
-          // Scale the base model viscosity value based on the porosity.
-          for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
+          const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");


Not something to modify in this PR, but just a self-reminder that this line of code (and others noted below) will apply to most of the reaction models and we might want to think about ways to avoid duplication in future PRs.

naliboff · 2024-06-06T21:16:58Z

source/material_model/reactive_fluid_transport.cc

+              for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
+                {
+                  const double porosity = std::max(in.composition[q][porosity_idx],0.0);
+                  out.viscosities[q] *= (1.0 - porosity) * exp(- alpha_phi * porosity);


Another line of code we might want to think about generalizing for multiple reaction models in future PRs (again, not a change for this PR).

@naliboff This is actually slightly different from what the Katz2003 model does. If everything was the same, it would be super easy to combine. But I agree, it would be nice not having to duplicate this.

naliboff · 2024-06-07T03:00:06Z

source/material_model/reactive_fluid_transport.cc

-                if (in.composition[q][porosity_idx] + porosity_change < 0)
-                  porosity_change = -in.composition[q][porosity_idx];
+          // Fill reaction rate outputs if the model uses operator splitting.
+          // Specifically, change the porosity (representing the amount of free water)


Suggested change

// Specifically, change the porosity (representing the amount of free water)

// Specifically, change the porosity (representing the amount of free fluid)

naliboff · 2024-06-07T03:00:23Z

source/material_model/reactive_fluid_transport.cc

-                  porosity_change = -in.composition[q][porosity_idx];
+          // Fill reaction rate outputs if the model uses operator splitting.
+          // Specifically, change the porosity (representing the amount of free water)
+          // based on the water solubility and the water content.


Suggested change

// based on the water solubility and the water content.

// based on the fluid solubility and the fluid content.

This and the above comment are regarding code not introduced in this PR, but I think it would be worth changing it now if I indeed this part of the code is not specific to one of the water-based reaction models (it looks like this is not the case).

jdannberg

Looks good to me!
The only test that's not working is the one you added. So if you've made sure any differences between the rising melt blob freezing and viscoplastic melt blob freezing results are due to the different solid viscosity, you can just updated the test case.

jdannberg · 2024-06-08T18:08:44Z

source/material_model/reactive_fluid_transport.cc

+              for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
+                {
+                  const double porosity = std::max(in.composition[q][porosity_idx],0.0);
+                  out.viscosities[q] *= (1.0 - porosity) * exp(- alpha_phi * porosity);


@naliboff This is actually slightly different from what the Katz2003 model does. If everything was the same, it would be super easy to combine. But I agree, it would be nice not having to duplicate this.

naliboff

@Grant-Block - This is ready to go from my end as well, aside from one minor comment change in the test prm file and double checking that the change in test output is only due to the change in solid viscosity.

naliboff · 2024-06-08T20:35:21Z

tests/viscoplastic_melt_blob_freezing.prm

+#########################################################
+# This is a test for the melting rate functionality,
+# in particular the freezing of melt.
+# It uses the melt simple model.


Suggested change

# It uses the melt simple model.

# It uses the reactive fluid transport material model composited

# on top of the visco plastic material model.

Hi @naliboff and @jdannberg, I double checked that the only difference in the output between this and the melt simple version of the test is the solid viscosity, so I applied the updated. Thank you!

naliboff · 2024-06-08T20:37:11Z

/rebuild

jdannberg reviewed Jun 5, 2024

View reviewed changes

Added melt outputs section to reactive fluid transport for katz2003 m…

61d1e34

…odel

Grant-Block force-pushed the add-katz2003-rft branch from 4ef9186 to 61d1e34 Compare June 6, 2024 20:01

jdannberg reviewed Jun 6, 2024

View reviewed changes

naliboff reviewed Jun 7, 2024

View reviewed changes

Grant-Block added 4 commits June 6, 2024 21:28

Improved code documentation and organization. Added a changelog file.

36fa765

Deleted header to try to fix spellcheck

d01aa87

Added header back to screen-output

b6fb9bb

Fixed test output

38c6c64

jdannberg approved these changes Jun 8, 2024

View reviewed changes

jdannberg added the ready to test label Jun 8, 2024

naliboff approved these changes Jun 8, 2024

View reviewed changes

Minor prm documentation change and updated test results

f6f66b9

jdannberg merged commit 122a0b2 into geodynamics:main Jun 9, 2024
8 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Added reaction and melt output calculations to reactive fluid flow model #5795

Added reaction and melt output calculations to reactive fluid flow model #5795

Grant-Block commented Jun 5, 2024

Grant-Block commented Jun 5, 2024

jdannberg commented Jun 5, 2024

jdannberg left a comment

jdannberg Jun 5, 2024

Grant-Block Jun 5, 2024

jdannberg Jun 5, 2024

Grant-Block commented Jun 6, 2024

jdannberg left a comment

jdannberg Jun 6, 2024

jdannberg Jun 6, 2024

jdannberg Jun 6, 2024

Grant-Block commented Jun 7, 2024

naliboff left a comment

naliboff Jun 6, 2024

naliboff Jun 6, 2024

jdannberg Jun 8, 2024

naliboff Jun 7, 2024

naliboff Jun 7, 2024

naliboff Jun 7, 2024

jdannberg left a comment

jdannberg Jun 8, 2024

naliboff left a comment

naliboff Jun 8, 2024

Grant-Block Jun 8, 2024

naliboff commented Jun 8, 2024

		"The reference temperature $T_0$. The reference temperature is used "
		"in both the density and viscosity formulas. Units: \\si{\\kelvin}.");

	// Specifically, change the porosity (representing the amount of free water)
	// Specifically, change the porosity (representing the amount of free fluid)

	// based on the water solubility and the water content.
	// based on the fluid solubility and the fluid content.

	# It uses the melt simple model.
	# It uses the reactive fluid transport material model composited
	# on top of the visco plastic material model.

Added reaction and melt output calculations to reactive fluid flow model #5795

Added reaction and melt output calculations to reactive fluid flow model #5795

Conversation

Grant-Block commented Jun 5, 2024

Before your first pull request:

For all pull requests:

For new features/models or changes of existing features:

Grant-Block commented Jun 5, 2024

jdannberg commented Jun 5, 2024

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

Grant-Block commented Jun 6, 2024

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

Grant-Block commented Jun 7, 2024

naliboff 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

jdannberg left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

naliboff commented Jun 8, 2024