(WIP) Add initial implementation of Maxwell-Wiechert rheology with two Maxwell arms by alarshi · Pull Request #7228 · geodynamics/aspect

alarshi · 2026-07-30T07:48:40Z

This PR is largely based on elastic rheology but with modification to the rhs terms following solution to several parallel maxwell arms (Maxwell wiechert elastic implementation by @bobmyhill )
@bobmyhill : This is much more simplified than what you derived, but it would be helpful to have your thoughts on the computation (ln 681–721 in the .cc file), everything else is related to prms that are needed for the computation of this term.

Before your first pull request:

I have read the guidelines in our CONTRIBUTING.md document.

For all pull requests:

I have followed the instructions for indenting my code.
I followed the AI Policy:
- Significant parts of this PR were written by AI (please add a sentence below).
- AI has not been used in significant ways.

If yes, please describe your usage of AI models in the creation of this pull request

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.

Add changes to stress compositional fields. Add stress fraction as compositional field. Restructure the reaction outputs.

bobmyhill

Hi @alarshi ,

Thanks for making a start on this cool project!

This PR is currently quite difficult to review, because a lot is derived from the elasticity rheology. Could you start with a commit that duplicates elasticity.cc (and .h) into maxwell_wiechert.cc (and .h) and just change the class name and initialisation? Then the github review system will be easier to use and you'll be able to track which things you changed in each commit.

After that, you could add what you've done here on top of that commit (if you are happy that you can follow the changes you made so far), or make smaller commits, slowly adding up to the model that you want to have. For example:

Commit 2: add the "stress fraction" fields (n - 1 fields for n Maxwell arms) and read them in to the rheology model. In the first commit you could have these fields do nothing and just be advected passively. You don't need to add a "Stress fractions" parameter in the rheology model - we will just need to strip these out again later.
Commit 3: Add a test for the stress fraction outputs to make sure that they work as you intend.
Commit 4: Add the stress fraction reaction terms. Again, they can just be placeholders in this commit - you'll have a later commit that adds their real values.
Commit 5: ...

I also gave you access to the overleaf that we're using to build this rheology model; I think that when we've populated the table with all the equivalencies between elasticity and maxwell_wiechert it will be a lot easier for you to substitute in the changes required to make this rheology work.

bobmyhill · 2026-08-04T14:03:15Z

+  <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _aspect_material_model_rheology_maxwell_wierchert_h


Suggested change

#ifndef _aspect_material_model_rheology_maxwell_wierchert_h

#ifndef _aspect_material_model_rheology_maxwell_wiechert_h

bobmyhill · 2026-08-04T14:03:23Z

+*/
+
+#ifndef _aspect_material_model_rheology_maxwell_wierchert_h
+#define _aspect_material_model_rheology_maxwell_wierchert_h


Suggested change

#define _aspect_material_model_rheology_maxwell_wierchert_h

#define _aspect_material_model_rheology_maxwell_wiechert_h

bobmyhill · 2026-08-04T14:04:22Z

+      };
+
+      template <int dim>
+      class MaxwellWierchert : public ::aspect::SimulatorAccess<dim>


Suggested change

class MaxwellWierchert : public ::aspect::SimulatorAccess<dim>

class MaxwellWiechert : public ::aspect::SimulatorAccess<dim>

bobmyhill · 2026-08-04T14:05:41Z

      else if (input == "entropy")
        return CompositionalFieldDescription::entropy;
+      else if (input == "stress fraction")
+        return CompositionalFieldDescription::entropy;


Suggested change

return CompositionalFieldDescription::entropy;

return CompositionalFieldDescription::stress_fraction;

bobmyhill · 2026-08-04T14:06:53Z

+      {
+        prm.declare_entry ("Elastic shear moduli", "75.0e9",
+                           Patterns::List(Patterns::Double(0.)),
+                           "List of elastic shear moduli, $G$, for each of the maxwell arms. "


Suggested change

"List of elastic shear moduli, $G$, for each of the maxwell arms. "

"List of elastic shear moduli, $G$, for each of the Maxwell arms. "

bobmyhill · 2026-08-04T14:13:27Z

+                           "is equal to the elastic time step. The default value of this parameter is "
+                           "'unspecified', which throws an exception during runtime. In order for "
+                           "the model to run the user must select 'true' or 'false'.");
+        prm.declare_entry ("Fixed elastic time step", "1.e3",


bobmyhill · 2026-08-04T14:13:50Z

+                           "timestep; afterwards on if 'Used fixed elastic time step' is true. "
+                           "Units: \\si{\\year} if the 'Use years instead of seconds' parameter is set; "
+                           "\\si{\\second} otherwise.");
+        prm.declare_entry ("Stabilization time scale factor", "1.",


bobmyhill · 2026-08-04T14:14:12Z

+      {
+
+        AssertThrow(this->get_parameters().enable_elasticity == true,
+                    ExcMessage("Rheology model maxwell wierchert only works if 'Enable elasticity' is set to true"));


Spelling and capitalization here and elsewhere.

bobmyhill · 2026-08-04T15:30:48Z

+
+                const double dtc = calculate_timestep_ratio() * elastic_timestep();
+
+                // Everything bellow and above this section is modified from the original rheology of Elasticiy


below and above? so, everything is modified?

bobmyhill · 2026-08-04T15:33:24Z

+
+                for (unsigned int a = 0; a < number_of_maxwell_arms; ++a)
+                  {
+                    const unsigned int f_index = this->introspection().compositional_index_for_name("stress_fraction_" + Utilities::int_to_string(a+1));


So the fractions are fixed? I don't think that this will do what you want it to do. I suspect that this model is equivalent to a single arm.

alarshi added 2 commits July 30, 2026 08:54

Add initial implementation with two maxwell arms.

6fad2a9

Add changes to stress compositional fields. Add stress fraction as compositional field. Restructure the reaction outputs.

Remove files that will be another PR.

9e94efd

bobmyhill changed the title ~~(WIP) Add initial implementation with two maxwell arms.~~ (WIP) Add initial implementation of Maxwell-Wiechert rheology with two Maxwell arms Jul 30, 2026

bobmyhill reviewed Aug 5, 2026

View reviewed changes

Navigation Menu

Uh oh!

(WIP) Add initial implementation of Maxwell-Wiechert rheology with two Maxwell arms - #7228

(WIP) Add initial implementation of Maxwell-Wiechert rheology with two Maxwell arms#7228
alarshi wants to merge 2 commits into
geodynamics:mainfrom
alarshi:test_simple_implementation

alarshi commented Jul 30, 2026

Uh oh!

bobmyhill left a comment

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

bobmyhill Aug 4, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

	#ifndef _aspect_material_model_rheology_maxwell_wierchert_h
	#ifndef _aspect_material_model_rheology_maxwell_wiechert_h

	#define _aspect_material_model_rheology_maxwell_wierchert_h
	#define _aspect_material_model_rheology_maxwell_wiechert_h

	class MaxwellWierchert : public ::aspect::SimulatorAccess<dim>
	class MaxwellWiechert : public ::aspect::SimulatorAccess<dim>

	return CompositionalFieldDescription::entropy;
	return CompositionalFieldDescription::stress_fraction;

	"List of elastic shear moduli, $G$, for each of the maxwell arms. "
	"List of elastic shear moduli, $G$, for each of the Maxwell arms. "


		const double dtc = calculate_timestep_ratio() * elastic_timestep();

		// Everything bellow and above this section is modified from the original rheology of Elasticiy

Uh oh!

Conversation

alarshi commented Jul 30, 2026

Before your first pull request:

For all pull requests:

For new features/models or changes of existing features:

Uh oh!

bobmyhill left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants