I would like to use the pylith_powerlaw_gendb.py function to create a spatial database for the power law rheology. In the doc, it says that we can make all the information in a file called powerlaw_gendb.cfg and that we need input spatial databases with the laboratory parameters. But I haven’t found any examples of how to use this function and I have no idea what these files look like.

We do not have an example parameter file for the utility pylith_powerlaw_gendb.py in the PyLith code base. This has been on our todo list and we haven’t gotten around to it. It looks like the utility needs a little updating to work with the power-law bulk rheology implementation in PyLith v3. I will try to create a working example in the next week.

Thank you very much to have created an example so quickly ! I have dowloaded the different files needed but when I try the utility function pylith_powerlaw_gendb.py as it is written in the documentation (without changing anything just to check), I get this error message :

I have tried with the SimpleDB spatial parameter, but I still have the same error …
Do I need to activate a specific thing to use the ‘pylith_powerlaw_gendb.py’ function ?

Thank you for your help, I have managed to run the simulation with an isotropic powerlaw rheology.
However, I have a question is there any threshold for parameters in the powerlaw rheology ?
For parameters log_flow_constant activation_energy (kJ/mol) power_law_exponent flow_constant_scale such like : 3.7 100 4.5 6.0 it is running but if I change the activation energy to 52 (laboratory observation for rocksalt) I have an error at the time step 2.

The error is
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Work array still checked out
[0]PETSC ERROR: WARNING! There are option(s) set that were not used! Could be the program crashed before they were used or a spelling mistake, etc!
[0]PETSC ERROR: Option left: name:-fieldsplit_displacement_ksp_type value: preonly
[0]PETSC ERROR: Option left: name:-fieldsplit_displacement_pc_type value: lu
[0]PETSC ERROR: Option left: name:-fieldsplit_lagrange_multiplier_fault_ksp_type value: preonly
[0]PETSC ERROR: Option left: name:-fieldsplit_lagrange_multiplier_fault_pc_type value: lu
[0]PETSC ERROR: Option left: name:-ksp_converged_reason (no value)
[0]PETSC ERROR: See FAQ — PETSc 3.19.5 documentation for trouble shooting.
[0]PETSC ERROR: Petsc Development GIT revision: v3.18.0-56-g437e83fba6f GIT Date: 2022-10-12 17:01:26 -0400
[0]PETSC ERROR: /data/cycle/meyerpa/pylith/pylith-3.0.3-linux-x86_64/bin/mpinemesis on a named ist-calcul23 by meyerpa Wed Sep 27 10:26:29 2023
[0]PETSC ERROR: Configure options --prefix=/opt/pylith/dist --with-c2html=0 --with-x=0 --with-clanguage=C --with-mpicompilers=1 --with-shared-libraries=1 --with-64-bit-points=1 --with-large-file-io=1 --download-chaco=1 --download-f2cblaslapack=1 --download-ml --with-fc=0 --with-hwloc=0 --with-ssl=0 --with-x=0 --with-c2html=0 --with-lgrind=0 --with-hdf5=1 --with-zlib=1 --LIBS=-lz --with-debugging=0 --with-fc=0 CFLAGS=“-g -O2” CXXFLAGS=“-g -O2 -DMPICH_IGNORE_CXX_SEEK” FCFLAGS= CPPFLAGS=“-I/opt/pylith/dist/include -I/opt/pylith/dist/include” LDFLAGS=“-L/opt/pylith/dist/lib -L/opt/pylith/dist/lib64 -L/opt/pylith/dist/lib” PETSC_DIR=/opt/pylith/build/cig/petsc-pylith PETSC_ARCH=arch-pylith
[0]PETSC ERROR: #1 DMDestroy() at /opt/pylith/build/cig/petsc-pylith/src/dm/interface/dm.c:718
[0]PETSC ERROR: #2 TSDestroy() at /opt/pylith/build/cig/petsc-pylith/src/ts/interface/ts.c:2720
[0]PETSC ERROR: #3 void pylith::problems::TimeDependent::deallocate()() at …/…/…/pylith-3.0.3/libsrc/pylith/problems/TimeDependent.cc:92
Is there anything I can do to solve it ?

I suspect the error about “Object is in wrong state” is a secondary error and the real error is that the solver did not converge. This could be due to too large of a time step or insufficient spatial resolution (too coarse of a mesh). There are not built-in thresholds for powerlaw rheology parameters. You must select an appropriate discretization size and solver parameters (such as the time step) to resolve the resulting deformation and allow the solver to converge.

Some strategies to resolve this problem include:

Select the appropriate basis order for your solution and auxiliary fields. The default basis order is 1 which means the field is represented as a function that varies linearly in space. A basis order of 2 represents the field as a quadratic function (polynomial includes terms with x, y, xy, x^2 and y^2).

Reduce the time step.

Examine the solution to see if the discretization size is small enough to resolve the deformation and refine the mesh appropriately.