Skip to content

Add Viscosity Limit to Reactive Fluid Transport#6883

Open
danieldouglas92 wants to merge 2 commits intogeodynamics:mainfrom
danieldouglas92:limit_viscosity_in_rft
Open

Add Viscosity Limit to Reactive Fluid Transport#6883
danieldouglas92 wants to merge 2 commits intogeodynamics:mainfrom
danieldouglas92:limit_viscosity_in_rft

Conversation

@danieldouglas92
Copy link
Contributor

In the evaluate function of Reactive fluid transport model, the base material model is first evaluated and then a fluid weakening term is multiplied to the viscosity computed in the base model. For base models like Visco plastic, this means that the viscosity that is returned by Reactive fluid transport can be much lower than the Minimum viscosity parameter set within the base model.

This PR adds a new parameter to the Reactive fluid transport model called Minimum weakened viscosity that sets a lower bound on the viscosity after applying this fluid weakening factor. I was hoping that there might be a better way to do this without defining an essentially duplicated parameter, but the problem is that not all base models contain a parameter that sets a minimum viscosity value (i.e. the Simple material model). If there is an alternative way to generalize this, I would love to hear other ideas!

@danieldouglas92
Copy link
Contributor Author

FYI @naliboff @xlia62

@danieldouglas92 danieldouglas92 force-pushed the limit_viscosity_in_rft branch 3 times, most recently from 34d16db to feb417c Compare March 6, 2026 03:48
Copy link
Member

@gassmoeller gassmoeller left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice, just two minor comments.

{
const double porosity = std::max(in.composition[q][porosity_idx],0.0);
out.viscosities[q] *= (1.0 - porosity) * std::exp(- alpha_phi * porosity);
const double weakend_viscosity = out.viscosities[q] * (1.0 - porosity) * std::exp(- alpha_phi * porosity);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
const double weakend_viscosity = out.viscosities[q] * (1.0 - porosity) * std::exp(- alpha_phi * porosity);
const double weakened_viscosity = out.viscosities[q] * (1.0 - porosity) * std::exp(- alpha_phi * porosity);

double shear_to_bulk_viscosity_ratio;
double min_compaction_visc;
double max_compaction_visc;
double min_weakened_visc;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know the other variable were around already, but can you change this to make it easier to read:

Suggested change
double min_weakened_visc;
double min_weakened_viscosity;

@danieldouglas92
Copy link
Contributor Author

@gassmoeller Thanks for the quick review and for catching that spelling mistake! I opened up a quick follow up PR to change the parameter name of the compaction viscosity limiters to be consistent with the new parameter in this pr

Copy link
Contributor

@naliboff naliboff left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@danieldouglas92 - Nice, thank you for putting this PR together! I only have to minor suggested changes on top of those already noted by @gassmoeller.

@@ -0,0 +1,32 @@
# This model tests the minimum viscosity within the reactive fluid transport model
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# This model tests the minimum viscosity within the reactive fluid transport model
# This model tests that the minimum viscosity within the reactive fluid transport model

# weakening factor of 30. Without a minimum value for the viscosity, this would
# result in the material viscosity of 4.481e17 Pa s, which is below the minimum
# value that is specified in the Visco Plastic base material model. By setting
# the parameter Minimum weakened viscosity to 1e10 Pa s, we ensure that the viscosity
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
# the parameter Minimum weakened viscosity to 1e10 Pa s, we ensure that the viscosity
# the parameter Minimum weakened viscosity to 1e19 Pa s, we ensure that the viscosity

@danieldouglas92 danieldouglas92 force-pushed the limit_viscosity_in_rft branch from e3ba727 to 7d33206 Compare March 6, 2026 22:16
@danieldouglas92
Copy link
Contributor Author

@naliboff Thanks for catching the typos! Just addressed your comments

@danieldouglas92 danieldouglas92 force-pushed the limit_viscosity_in_rft branch from 7d33206 to 0f45daf Compare March 8, 2026 17:28
@danieldouglas92
Copy link
Contributor Author

I'm not sure why the indent tester is failing? I've run make indent locally and no files are changed, and the indent artifact that the tester outputs is empty

@naliboff
Copy link
Contributor

naliboff commented Mar 9, 2026

/rebuild

@gassmoeller
Copy link
Member

The code looks all good now.

The tester is failing, because since #6881 we also check in this tester if the parameter documentation needs to be updated, and you added a new parameter here. You can update the parameter documentation manually by running ./update_parameters.sh ../build/aspect from the doc/ directory, or you can download a changes file from the tester. Unfortunately, the tester had a bug that didnt let it upload the changes file correctly. Can you take a look at #6902 and merge that (if it passes the tests)? Then rebase this PR and you should get a changes file prepared on the results page of the tester (just like for the indent and results tester).

@danieldouglas92
Copy link
Contributor Author

danieldouglas92 commented Mar 11, 2026

The first commit worked for models that I was running (which did not use the Katz2003 melting model), but @xlia62 pointed out that it wasnt working when using the katz2003 melting model. This is because the katz2003 melting model calculates this weakening factor on it's own within reaction_model/katz2003_melting_model.cc.

I pushed a second commit which removes this code from katz2003_melting_model.cc and moves the code in reactive_fluid_transport.cc around so that the weakening is applied regardless of the reaction model being used. I think this makes sense to avoid needing to potentially make individual viscosity limiters for each reaction model, and I think this should just work, but we'll see if the testers pass

@danieldouglas92
Copy link
Contributor Author

@gassmoeller Ah that makes sense! I'll push more changes in a minute with the updated parameter.

@danieldouglas92 danieldouglas92 force-pushed the limit_viscosity_in_rft branch 2 times, most recently from 20bbb07 to 614225d Compare March 11, 2026 16:31
@danieldouglas92 danieldouglas92 force-pushed the limit_viscosity_in_rft branch from 614225d to 8ea5508 Compare March 11, 2026 17:08
@danieldouglas92
Copy link
Contributor Author

danieldouglas92 commented Mar 11, 2026

The changes to katz2003 caused two of the tests (viscoplastic_melt_blob_freezing and shear_heating_with_melt) to fail. viscoplastic_melt_blob_freezing failing is confusing to me because it uses the reactive fluid transport model, but rising_melt_blob_freezing passes which uses the melt simple material mode. shear_heating_with_melt also uses the melt simple material model, and so it would make sense why this test might fail since the melt weakening is now only being calculated if the katz2003 reaction model uses the reactive fluid transport model.

Clearly this is why there was an if statement in reactive_fluid_transport that excluded the fluid weakening from being calculated if the reaction model was equal to katz2003. Because the katz2003 interfaced with material models outside of reactive fluid transport, the weakening had to be determined internally within the reaction model itself.

An easy solution to this would be to just add an identical parameter to the katz2003 reaction model that allows the user to limit the reduction of the shear viscosity, but this feels pretty hacky.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants