Issue with Compositional Fields when using Multicomponent Material Model with Ascii Initial Composition

ASPECT
1

Hi all!

I’m new to ASPECT and trying to recreate the initial conditions of this paper https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2015tc003908 (Flat-slab subduction, topography, and mantle dynamics in southwestern Mexico, Gérault et al. 2015), however I’ve run into a bit of an issue. To have the geometry shown in figure 2, I used the multicomponent Material model with an ascii data Initial composition model. There are 11 unique combinations of density and viscosity according to Figure 2, so I made an ascii file that assigned each point in the chunk to one of the compositional fields. (It is an instantaneous model with thermal expansivities and conductivities = 0)

However, after running ASPECT (version 3.1.0-pre) I noticed some strange things. In ParaView, I can see an Array for each of my 11 compositional fields, however only the first three are correctly assigned, the rest have range [0,0]. Additionally, when looking at the density and viscosity arrays, their ranges also only include the values from the first three compositional fields. Everything else, including the background mantle has not been incorporated.

image
image260×249 5.01 KB

Secondly, in the output that ASPECT prints to the command line, it skips all the compositional fields after the first three because “RHS is zero”:
*** Timestep 0: t=0 seconds, dt=0 seconds
Skipping temperature solve because RHS is zero.
Solving top_mantle system … 0 iterations.
Solving mid_mantle system … 0 iterations.
Solving bot_mantle system … 0 iterations.
Skipping subducting_plate composition solve because RHS is zero.
Skipping weak_zone composition solve because RHS is zero.
Skipping overriding_plate_1 composition solve because RHS is zero.
Skipping overriding_plate_2 composition solve because RHS is zero.
Skipping overriding_crust composition solve because RHS is zero.
Skipping subducting_crust composition solve because RHS is zero.
Skipping weak_line_1 composition solve because RHS is zero.
Skipping weak_line_2 composition solve because RHS is zero.
Solving Stokes system (GMG)… 20+0 iterations.
Postprocessing:
Writing graphical output: output-flat_slab/solution/solution-00000
RMS, max velocity: 4.61e-17 m/s, 1.59e-16 m/s
Termination requested by criterion: end time

At first I thought there might have been a limit to only three compositional fields, however I swapped two around and still only the original first three were incorporated:
Skipping temperature solve because RHS is zero.
Solving top_mantle system … 0 iterations.
Solving mid_mantle system … 0 iterations.
Skipping subducting_plate composition solve because RHS is zero.
Solving bot_mantle system … 0 iterations.
…

I’ve checked the ascii file and there are definitely points that are assigned to the other compositions. I also tried increasing the global refinement level to 7 in case it was previously too coarse to recover the slabs, however it didn’t make a difference. I also ran it for 200,000 years instead of instantaneously but the same issue persisted.

Has anyone encountered this before, or have any advice on how to fix it?

Thank you!

2

@acol-63 - Welcome and thank you for posting your question to the forum.

I think the issue likely lies within the script used to generate the ascii data, but a minor adjustment to the PRM could also be needed. Can you provide both so we can help you debug?

As an aside, for these types of complex geometries I highly recommend using the Geodynamic World Builder. It will make generating complex geometries easier and allow you to capture key features at much higher resolution as it provides a mesh-independent representation that is re-queried with each successive initial adaptive refinement step.

Cheers,

John

3

Hi John,

I made the ascii file by:

  1. Using WebPlotDigitiser to extract the boundaries of each region into .txt files of points
  2. Ran combine_regions.sh to combine all the polygons into one file
  3. In GMT, used grdmask and then grd2xyz to turn the polygons into an ascii grid (annulus_grid.gmt)
  4. Realised the ascii file had to be in a specific format so I used columns.sh to change it and check the formatting

Then I put the output of that into my .prm file (gerault.prm).

The scripts are here:

Scripts.zip (3.1 KB)

Thank you for your reccomendation, I will take a look at the world builder.

Thanks for your reply!

4

Hi @acol-63,

Thank you for sending the files over. I looked at the PRM and did not see any obvious issues.

A good first check would be to pick a point where you know an existing compositional field should be present (value = 1) but has a value of 0, and then backtrack through each step to see where the things are going wrong.

Cheers,

John

5

Hi @acol-63 ,

Welcome to ASPECT. I have a simple way to figure out which and where goes wrong: Adding compositional field one by one and see what’s happening, so that you can pinpoint where is problem is.

Cheers,
Mingming

6

Hi John,

I’ve checked several points in the ascii file and they all appear to have 1 in the correct column for that point and 0 in the rest. Would it be helpful if I attached the ascii file here?

Thanks

7

Hi Mingming,

I removed all the fields apart from one that was not one of the 3 that seem to be working (I chose subducting_plate). I checked in the ascii file and there are 563259 points (~1.56% of them) that should have the subducting_plate composition.

However, when I ran ASPECT I once again received the “Skipping subducting_plate composition solve because RHS is zero” line in the terminal output. When I checked in ParaView, the range for that composition was zero.

image

And the viscosity and density everywhere were the reference values I set in the .prm file.

I’m not really sure what could be causing it to ignore the ascii file for the 8 compositions. Do you have any recommendations?

Thanks for your reply.

8

Hi @jbnaliboff,

Here is the ASCII file: https://drive.google.com/file/d/1yEvxOp6Qh352lDX8Ba0sx-jIht1M3E5o/view?usp=sharing

Sorry it was too large to attatch here so I’ve sent it as a Google Drive link.

Thanks

9

Hi @acol-63,

Thank you for sharing your input files! Just to add to what @jbnaliboff and @optimux suggested, the reason you’re missing the other compositional fields is because ASPECT expects the ascii grid in spherical coordinates to have values in the order r (meters), phi (radians), theta (radians), followed by the compositional field values. However, your input file specifies the second column, i.e., phi, in degrees. Therefore, some compositional fields did not show up even though they were present because they were interpreted as radians.

I was going to point you to the relevant section in the manual, but it turns out we don’t actually specify the units there. I will open an issue and address it in a subsequent pull request.

Regards,
Arushi

10

Hi Arushi,

I changed my phi column to radians and it’s working now! I thought that since the chunk geometry min/max longitude is in degrees, the ascii file should also be in degrees so I didn’t think to change it, even though I probably should have.

Thanks for your help!

11

Hi @acol-63 ,

So sorry for this late reply, I was too busy then. It looks like @arushi_saxena solved your problem, that’s good. If you have more questions, I’m happy to help.

Cheers,
Mingming