Input file issue for the Continental Rift model

zch_aspect_200km_no_fastscape.prm (21.6 KB)
log.txt (182.7 KB)
Dear Experts,

I am currently using the Ubuntu 23.1 system and have recently downloaded the ASPECT 2.6-pre version code for numerical simulation of continental rifting. The versions of other software I am using are deal.II 9.5.2, Trilinos 13.2.0, and p4est 2.3.2.

I have created an inputfile based on the prompts from ASPECT, which can be recognized and iterated. However, the solver quickly fails to converge. I have tried adjusting many parameters, such as rheological parameters, temperature field parameters, and solver parameters, but none have resolved the convergence issue. I noticed that in the output files, velocity and strain rates are concentrated at the boundary between the lithosphere and asthenosphere, which seems incorrect, but I am unsure how to modify it. Therefore, I have uploaded my input parameters and the log.txt file. Could you please advise me on how to resolve this issue?

Additionally, I found that the continental_extension.prm in cookbooks can only use the single Advection, iterated Newton Stokes solver, while the single Advection, iterated defect correction Stokes and single Advection, iterated Stokes solvers cannot iterate smoothly. Is this phenomenon common?

Thank you for your assistance.

Best regards,

Chenghui Zeng

Hi Chengzui,

Thank you for posting to the forum and providing a PRM file. Some initial replies below:

  1. For the first issue you described, did you start with the continental extension cookbook and then modify to the current setup? If yes, do you recall at what stages you started have convergence issues? Going back to a setup that works, and then systematically introducing changes to see where the error arises may the easiest path forward. However, an initial suggestion below

  2. Again in regards to the first issue, it looks like you are using a different flow law between the asthenosphere and lithospheric mantle? What is the type of flow law you are using for the asthenosphere? If these flows produce a large viscosity contrasts across the LAB, then it is not so surprising that you would get very high strain rates there. Your minimum viscosity is also sufficiently low (1e18 Pa s) that you may hit that minimum value if the asthenosphere flow law is quite weak. In practice, often times my colleagues and I have often had the majority of convergence issues in 2D/3D rift (or subduction/collision) models arise from sharp viscosity gradients in the asthenosphere. Using smaller time steps can help, but in many cases we had alter the asthenospehere viscous flow laws and/or increase the minimum viscosity. These are at least a few ideas for things to test.

I found that the continental_extension.prm in cookbooks can only use the single Advection, iterated Newton Stokes solver, while the single Advection, iterated defect correction Stokes and single Advection, iterated Stokes solvers cannot iterate smoothly. Is this phenomenon common?

That is new, at least to my knowledge. Convergence in the past is generally very stable (converges in 2-5 iterations) after the first time step. Can you make a plot of the nonlinear convergence behavior you observed between these three settings? I actually have a new version of the continental extension cookbook nearly ready to submit, which uses the Newton solver and adds elasticity. If it is helpful, here are the solver settings I have been using there:

set Nonlinear solver scheme                = single Advection, iterated Newton Stokes
set Nonlinear solver tolerance             = 1e-5


subsection Solver parameters
  subsection Stokes solver parameters
    set Stokes solver type = block AMG
    set Number of cheap Stokes solver steps = 2000
    set Linear solver tolerance                        = 1e-7
    set GMRES solver restart length               = 100
    set Use full A block as preconditioner       = true
  end

  subsection Newton solver parameters
    set Max pre-Newton nonlinear iterations          = 10
    set SPD safety factor                                       = 0.9
    set Nonlinear Newton solver switch tolerance = 1e-7
    set Max Newton line search iterations             = 5
    set Maximum linear Stokes solver tolerance   = 1e-5
    set Use Newton residual scaling method       = false
    set Use Newton failsafe                                 = true
    set Stabilization preconditioner                     = SPD
    set Stabilization velocity block                       = SPD
    set Use Eisenstat Walker method for Picard iterations = false
  end
end

Please be sure that you are a version of ASPECT compiled from the main development branch, as a number of improvements to the Newton solver were recently introduced during the summer.

I hope this is helpful.

Cheers,
John

Dear John,

I hope this message finds you well.

First and foremost, I would like to express my sincere gratitude for your detailed and professional response. Your suggestions have been incredibly helpful in resolving my initial issues. By adjusting the minimum viscosity limit to 1e20, my model has been able to iterate smoothly. The significant viscosity gradient between the asthenosphere and lithospheric mantle, due to the use of wet olivine material at a temperature of 1380°C, was indeed a critical factor. I will continue to conduct rheological tests to refine this aspect of my model.

Additionally, the new solver parameters you provided have notably accelerated the iteration rate of my model. The issue with the cookbooks has been resolved, and I am truly appreciative of your assistance.

As I am still in the process of learning ASPECT and familiarizing myself with its various functionalities, I have a few more questions that I would like to seek your guidance on:

  1. Displaying Background Field in ParaView:
  • In ParaView, I am able to visualize compositional fields such as upper_crust, lower_crust, and mantle_lithosphere. However, I am unsure how to output the background field. Could you advise on which parameters to include in the postprocess section to visualize the background field? I need to display the background information in ParaView for creating a state.
  1. Using Multiple Material Model Plugins:
  • Within the material_model, there are various plugins available. Is it possible to use these plugins simultaneously? For instance, I am interested in calculating the melt fraction of the background (asthenosphere) during continental rifting. I am considering using plugins like latent heat melt or melt simple, and then visualizing the melt fraction in the postprocess section with visualization. However, there is no option for List model of names in the material_model. Does this imply that it is not possible to enable multiple plugins at once?
  1. Melt Transport and Solver Configuration:
  • When using melt functionalities, it is common to enable the melt transport feature, which requires setting set Use operator splitting = true. However, this seems to prevent the use of the Newton solver. Could you provide guidance on alternative solvers, such as single Advection, iterated Stokes and iterated Advection and Stokes, and any considerations I should keep in mind when selecting these solvers?

I eagerly await your expert advice on these matters. Your assistance is invaluable, and I greatly appreciate your time and effort.

Thank you once again,

Best regards,
Chenghui Zeng

Hi @zengchenghui,

My pleasure and I’m glad the previous suggestions were helpful.

Answers to the additional questions below.

Could you advise on which parameters to include in the postprocess section to visualize the background field?

There is actually no field or information that is currently output specifically for the background field.

In fact, it is not a field that is advected, but rather a field that is used to account for cases when the other fields representing lithologies do not sum up to a volume fraction of 1 at a given point.

In this case, the background field takes up the remaining volume fraction and its assigned properties are used accordingly when the properties between fields are averaged. If the background field has different values than the other fields, this should be reflected in the output material properties.

You can also obtain the relative volume fraction of the background field through post processing by summing the values of the compositional fields representing lithologies at every point, and then subtracting this value from 1. The remainder will represent the background field.

However, there is no option for List model of namesin thematerial_model . Does this imply that it is not possible to enable multiple plugins at once?

See this test case for an example of how material model compositing can work.

Also, here is a link to the compositing material model description and the relevant section in the manual.

When using melt functionalities, it is common to enable the melt transport feature, which requires setting set Use operator splitting = true . However, this seems to prevent the use of the Newton solver.

Yes, at present the melt functionality is not integrated with the Newton solver. You will need to use indeed need to use single Advection, iterated Stokes and iterated Advection and Stokes. The general consideration to keep in mind is to make sure the linear and nonlinear solvers are sufficiently strict to achieve a reasonable level of accuracy for the problem at hand.

Cheers,
John

Dear Experts,

I hope this email finds you well.

I am currently using the Ubuntu 23.1 system and have recently downloaded the ASPECT 2.6-pre version code for numerical simulation of continental rifting. The versions of other software I am using are deal.II 9.5.2, Trilinos 13.2.0, and p4est 2.3.2.

I would like to apologize for deleting my previous inquiry, as I wasn’t fully clear on the logic of the code and the questions I needed to ask. After some time spent testing, I have clarified the issues I am encountering. Below are my questions:

Integration of Melt Functions: I have merged functions from melt_simple.cc and katz2003_mantle_melting.cc into visco_plastic.cc so that I can simultaneously calculate the lithosphere’s rheological processes and changes in melting degree. This integration allows for easier incorporation of melting effects on density and viscosity in the future rather than in postprocess stage. Currently, I have simply combined all functions, and the output parameters such as density and viscosity (out.densities, out.viscosities) are not yet affected by the melting degree. If possible, could you review the evaluate function in my code and provide suggestions on how to relate the melting degree to density and viscosity?

Background Field Output: I am trying to output a background field to visualize it in Paraview software. I added an asthenosphere field in the output function and assigned it the value of volume_fractions[0] (as noted in the commented code in the first image). Although the compilation of ASPECT is successful, I encounter an error during the iteration (as shown in the second image). It seems that the array I created is not initialized properly. If convenient, could you check my visco_plastic_grain.cc file and advise me on how to modify the assignment process for this array to allow smooth iteration and output of the background field?

Fastscape Plugin Parameters: To simulate the rift sedimentation process, I have coupled the Fastscape plugin. However, as shown in images three and four, even though plastic strain accumulates, there is no significant fluctuation in the topography, and the sediment thickness at 11 Ma appears abnormal, with faults not controlling the sedimentation process. This might be due to insufficient settings for the Fastscape parameters in my .prm file. Could you please guide me on how to set the Fastscape parameters correctly to enable ASPECT to simulate the sediment filling process accurately?

Melt Transport Functionality: When integrating melt_simple.cc with visco_plastic.cc, I retained the melt transport functionality. I believe this is a valuable feature; however, enabling it leads to a noticeable decrease in iteration speed. In the context of considering melting effects during continental rifting processes, should I enable this functionality? Are there any detailed references discussing the impacts of melt transport functionality? I have only found information in the manual concerning the effects of melt transport on large-scale tectonics, such as mantle plumes.

Lastly, I would like to mention that I have provided the modified plugin polygon_temp_comp_topo.zip, which is based on the work of Anne Glerum, Derek Neuharth and other experts. I have changed the plugin’s name, added middle crust compositions, and modified the temperature field calculation method.

I look forward to your professional guidance on my questions. Thank you very much for your assistance!

Best regards,

Chenghui Zeng





polygon_temp_comp_topo.zip (6.3 MB)
visco_plastic_melt.zip (6.2 MB)

20240823.prm (31.8 KB)

@zengchenghui Q&A forums like this do not work well with questions that have multiple parts. Would you mind opening separate questions for each part you are asking about?

Best
W.

@zengchenghui - Following on @bangerth comment, I would also encourage you to take a look at the following merged pull requests that began the process of enabling combining different solid material models with existing melt functionality. These may simplify your workflows and reduce the amount of new code you need to develop.

I’m sorry that I don’t understand this form of question and answer, and thank you for your suggestion. I will open a new question page and ask about the problems I have encountered one by one.

Best regards,

Chenghui Zeng

Thank you very much for your suggestions. I will carefully read the materials you provided.

Best regards,

Chenghui Zeng