contribute cookbook: simple subduction initiation #5271

cedrict · 2023-07-12T01:29:36Z

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.

Describe what you did in this PR and why you did it.

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:

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 for adding this!

I have two general comments in addition to the ones I left in the code:

(1) Can you add a test case? It can just include your input file (add the line include $ASPECT_SOURCE_DIR/cookbooks/subduction_initiation/subduction_initiation_compositional_fields.prm) and then maybe reduce the resolution/end time. This way, we make sure it will continue to work and the results will still be correct in the future. I would actually do this first (before any of the other changes I suggested) because then you can use it as well to make sure that changing the input file does not change your solution.

(2) This is more a philosophical question: Your water viscosity is 1e19, which is the same as the lowest value of the asthenosphere viscosity (case 6). Do you think that is still a reasonable approximation? Could you maybe add a line to the text on that?

Thank you!

jdannberg · 2023-07-12T23:21:24Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+lithosphere ($\rho_l=3300~\text{ kg m}^{-3}$, $\eta_l$) and
+asthenosphere ($\rho_a=3200~\text{ kg m}^{-3}$, $\eta_a$).
+Note that although the water viscosity is correct, this is probably not
+the value that was used in the authors since it would most likely lead to


either

Suggested change

the value that was used in the authors since it would most likely lead to

the value that was used by the authors since it would most likely lead to

or

Suggested change

the value that was used in the authors since it would most likely lead to

the value that was used in the original publication since it would most likely lead to

jdannberg · 2023-07-12T23:21:49Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+asthenosphere ($\rho_a=3200~\text{ kg m}^{-3}$, $\eta_a$).
+Note that although the water viscosity is correct, this is probably not
+the value that was used in the authors since it would most likely lead to
+numerical errors. We then set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').


Suggested change

numerical errors. We then set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').

numerical errors. We here set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').

jdannberg · 2023-07-12T23:39:35Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+Two parameter files are provided for Case 1: one using
+compositional fields, and one using the particle-in-cell approach.


You can directly link to the files, I think it should look like this:

Suggested change

Two parameter files are provided for Case 1: one using

compositional fields, and one using the particle-in-cell approach.

Two parameter files are provided for Case 1: one

[using compositional fields](https://github.com/geodynamics/aspect/blob/main/cookbooks/subduction_initiation/subduction_initiation_compositional_fields.prm), and one

[using the particle-in-cell approach](https://github.com/geodynamics/aspect/blob/main/cookbooks/subduction_initiation/subduction_initiation_particle_in_cell.prm).

jdannberg · 2023-07-12T23:43:24Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+<img src="results3.*" width="95%" />
+
+Case 1 (obtained with the compositional fields approach) at
+(a) 0 Ma; (b) 12 Ma; (c) 21 Ma; (d) 31 Ma; (e) 41 Ma.


Can you add units? I am not going to ask you to redo all the figures (although having the units directly with the colorbar would be great!) but it should at least be in the caption.

added units in caption. also re-uploading figure bc original was too low resolution.

jdannberg · 2023-07-12T23:45:48Z

cookbooks/subduction_initiation/subduction_initiation_particle_in_cell.prm

+#subsection Solver parameters
+#  subsection Stokes solver parameters
+#    set Number of cheap Stokes solver steps = 0
+#  end
+#end


Can you remove this if it is not needed?

jdannberg · 2023-07-12T23:50:01Z

cookbooks/subduction_initiation/subduction_initiation_compositional_fields.prm

+subsection Initial composition model
+  set Model name = function
+  subsection Function
+    set Variable names      = x,y
+    set Function constants  = L0=300e3
+    set Function expression = if((x>=L0 && y<162e3),1,0) ;\
+                              if((x<=L0 && y>120e3 && y<=170e3),1,0) ;\
+                              if((x>=L0 && y>=162e3 && y<=172e3),1,0) ;\
+                              if((x<=L0 && y>170e3) || (x>=L0 && y>172e3),1,0)
+  end
+end


I get that you could change this by hand if you wanted to change the domain height to 140 km for cases 5 and 6. But you would have to change all the numbers by hand and it would be really easy to make an error. Could you add another function constant for the domain hight and then express all y values as the height minus depth? That way, the user would only have to change one number for the box height.

good idea, thx!

jdannberg · 2023-07-12T23:50:58Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+Two parameter files are provided for Case 1: one using
+compositional fields, and one using the particle-in-cell approach.
+Other cases can be run by changing the viscosities and/or the domain size accordingly in the
+parameter files.


Suggested change

parameter files.

parameter files. Note that changing the domain size does not just require changing the geometry model, but also the initial composition model.

jdannberg · 2023-07-12T23:53:34Z

cookbooks/subduction_initiation/subduction_initiation_particle_in_cell.prm

@@ -0,0 +1,104 @@
+#parameter file for replicating Matsumoto & Tomoda 1983


In order to not duplicate code, can you include the other input file and only list the parameters that are different, rather than copying the whole file (I think you only need the output directory, the section on compositional fields, and the postprocessing)? That way, if someone wants to change the setup to other cases, they only have one file to change, not both.

have created a 'subduction_initiation_basis.prm' file which contains the common prm code for both.

jdannberg · 2023-07-12T23:54:11Z

cookbooks/subduction_initiation/subduction_initiation_particle_in_cell.prm

+    set Data output format       = vtu
+    set List of particle properties = initial composition, velocity
+    set Particle generator name = random uniform
+    set Interpolation scheme = cell average #default


Suggested change

set Interpolation scheme = cell average #default

set Interpolation scheme = cell average

bangerth

I liked the "sticky water" analogy. Perhaps honey?

I had to look up the viscosity of honey, and found it to be on the order of 50 Pa.s. So your water with a viscosity of 1e19 Pa.s is really sticky!

One also finds this publication (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4745488/#:~:text=Viscosity%20of%20all%20the%20honey,54.0%20to%2089.0%20kJ%2Fmol.) wherein the authors also see whether irradiation with gamma rays changes the viscosity -- though I will admit that I cannot see why anyone would find themselves interested in the effects of radition on viscosity.

bangerth · 2023-07-13T00:24:09Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+
+The setup for this experiment originates in {cite:t}`matsumoto:tomoda:1983`.
+
+In this very early computational geodynamics paper the authors


Suggested change

In this very early computational geodynamics paper the authors

In this very early computational geodynamics paper, the authors

bangerth · 2023-07-13T00:24:54Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+at the gigantic fracture zone in the Northeastern
+Pacific by means of numerical simulation
+of contact between two kinds of viscous fluid of different density.
+They then solve the incompressible isothermal linear Stokes equations


Suggested change

They then solve the incompressible isothermal linear Stokes equations

They then formulate the incompressible isothermal linear Stokes equations

bangerth · 2023-07-13T00:25:32Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+lithosphere ($\rho_l=3300~\text{ kg m}^{-3}$, $\eta_l$) and
+asthenosphere ($\rho_a=3200~\text{ kg m}^{-3}$, $\eta_a$).
+Note that although the water viscosity is correct, this is probably not
+the value that was used in the authors since it would most likely lead to


Suggested change

the value that was used in the authors since it would most likely lead to

the value that was used in the paper since it would most likely lead to

already corrected (see Julianne's comment above)

bangerth · 2023-07-13T00:25:52Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+asthenosphere ($\rho_a=3200~\text{ kg m}^{-3}$, $\eta_a$).
+Note that although the water viscosity is correct, this is probably not
+the value that was used in the authors since it would most likely lead to
+numerical errors. We then set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').


Suggested change

numerical errors. We then set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').

numerical errors. We thus set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').

:-)

Honey?

sticky honey is a hard sell ... is there a planet covered in honey with plate tectonics?

bangerth · 2023-07-13T00:32:13Z

cookbooks/subduction_initiation/doc/subduction_initiation.md

+the value that was used in the authors since it would most likely lead to
+numerical errors. We then set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').
+Also, all viscosities in the paper
+are expressed in Poise with $1~\text{ Poise}=0.1~\text{ Pa s}$.


Ah, come on. Who thought that was a good idea?

bangerth · 2023-07-13T00:34:18Z

cookbooks/subduction_initiation/subduction_initiation_compositional_fields.prm

@@ -0,0 +1,98 @@
+#parameter file for replicating Matsumoto & Tomoda 1983


Suggested change

#parameter file for replicating Matsumoto & Tomoda 1983

# Parameter file for replicating Matsumoto & Tomoda 1983

cedrict · 2023-07-14T21:51:29Z

I have changed the water viscosity to 1e18. It does not change the results shown in the cookbook very much at all. But 1e18 is indeed safer than 1e19 in case 6. I have also added a link to Crameri paper.

cedrict · 2023-07-14T22:19:34Z

ready to review again. Still need to do a test

bangerth · 2023-07-14T22:31:04Z

This looks good. Do you want us to merge this now, or wait for the test?

cedrict · 2023-07-14T22:48:09Z

please merge, am now looking into making a test (could take a while, it's been a looong time :)

bangerth · 2023-07-14T22:53:29Z

There's a failing check:

--- a/cookbooks/subduction_initiation/subduction_initiation_basis.prm
+++ b/cookbooks/subduction_initiation/subduction_initiation_basis.prm
@@ -46,4 +46,3 @@ subsection Mesh refinement
   set Strategy                      = composition
   set Coarsening fraction           = 0
 end
-

Perhaps an empty line at the end of a file?

Other than that, good to merge!

contribute_cookbook

b24fd8a

cedrict force-pushed the add_subd_init_cookbook branch from 0579b4a to b24fd8a Compare July 12, 2023 03:08

tjhei added the ready to test label Jul 12, 2023

jdannberg reviewed Jul 13, 2023

View reviewed changes

bangerth reviewed Jul 13, 2023

View reviewed changes

addressed all comments

299eced

bangerth added the ready to merge label Jul 14, 2023

indent

fa466b3

naliboff mentioned this pull request Jul 15, 2023

add test subduction initiation cookbook #5329

Merged

5 tasks

gassmoeller approved these changes Jul 15, 2023

View reviewed changes

gassmoeller merged commit 4065456 into geodynamics:main Jul 15, 2023
6 checks passed

	the value that was used in the authors since it would most likely lead to
	the value that was used by the authors since it would most likely lead to

	the value that was used in the authors since it would most likely lead to
	the value that was used in the original publication since it would most likely lead to

	numerical errors. We then set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').
	numerical errors. We here set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').

		Two parameter files are provided for Case 1: one using
		compositional fields, and one using the particle-in-cell approach.

	parameter files.
	parameter files. Note that changing the domain size does not just require changing the geometry model, but also the initial composition model.

		@@ -0,0 +1,104 @@
		#parameter file for replicating Matsumoto & Tomoda 1983

	set Interpolation scheme = cell average #default
	set Interpolation scheme = cell average


		The setup for this experiment originates in {cite:t}`matsumoto:tomoda:1983`.

		In this very early computational geodynamics paper the authors

	They then solve the incompressible isothermal linear Stokes equations
	They then formulate the incompressible isothermal linear Stokes equations

	the value that was used in the authors since it would most likely lead to
	the value that was used in the paper since it would most likely lead to

	numerical errors. We then set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').
	numerical errors. We thus set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').

		@@ -0,0 +1,98 @@
		#parameter file for replicating Matsumoto & Tomoda 1983

	#parameter file for replicating Matsumoto & Tomoda 1983
	# Parameter file for replicating Matsumoto & Tomoda 1983

contribute cookbook: simple subduction initiation #5271

contribute cookbook: simple subduction initiation #5271

Conversation

cedrict commented Jul 12, 2023

Before your first pull request:

For all pull requests:

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

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

bangerth 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

cedrict commented Jul 14, 2023

cedrict commented Jul 14, 2023

bangerth commented Jul 14, 2023

cedrict commented Jul 14, 2023

bangerth commented Jul 14, 2023