Initialization of coupled flow and mechanics problems

[francoishamon]

I am opening this discussion to brainstorm about initialization for coupled flow and mechanics problems.

In the simulations that we have to run, we are often given an initial hydrostatic pressure field (and an initial composition) and we are asked to initialize the mechanical stresses such that the problem is at equilibrium when the simulation starts. Following discussions with Nicola and Josh, it seems the following initialization procedure would work well:

1- Compute or load the hydrostatic pressure field before the simulation starts

2- Disable inelasticity

3- Prior to the beginning of the simulation, solve an elastic problem that accounts for the initial pressure and composition.

In the context of the multiphase poromechanics solver, we can do that by solving the linear system:

delta x = - J^-1 R 

where the Jacobian looks like

     | J_uu J_up J_us |
J =  |        I       |
     |            I   |  

and where the residual has been modified appropriately to fix pressure and compositions at their initial values.

In GEOSX, this can be implemented as an initialization time step (time step size does not matter), and the modification of the Jacobian can be implemented with our field specifications functions. The linear solver may have to be modified to accomodate this Jacobian structure.

4- Enable inelasticity

5- Zero out the displacements computed during the initialization time step

6- Start the "real" simulation at time t = 0s

In terms of implementation:

• Step 1 is already in develop
• Steps 2, 4, and 5 are in PR Implemented a task for pre/post equilibration in coupled problems #1898
• Step 3 is implemented in a PR opened for discussion: PR Implemented initialization elastic time step for coupled problems #1927

I tested the procedure on realistic large problems and it seems to work well for realistic large cases.

Does this initialization procedure look ok for everyone?
I don't have a good vision of all the use cases, so I may be missing something (faults and fractures, etc)

@castelletto1 @joshua-white @rrsettgast @klevzoff @CusiniM @jhuang2601

[TotoGaz]

I do not understand if step 5 after step 4 implies the computation of any equilibrium.
Maybe it's a detail, maybe not.

However, it guess that steps 1 to 5 could be computed in another simulation, which outputs could be defined as input fields for the full big one.

I don't have a good vision of all the use cases, so I may be missing something (faults and fractures, etc)

So maybe not committing the solvers to one unique initialization recipe and letting the users develop their own mix is fine? 🧙 🤣

EDIT: Or maybe your SolidMechanicsStateReset can receive a Solvers node instead? Defining which fields you would like to be copied as initial conditions for the main solver.

[francoishamon]

However, it guess that steps 1 to 5 could be computed in another simulation, which outputs could be defined as input fields for the full big one.

Yes, definitely, there is always the possibility to output a restart file at the end of initialization, and then run the "real" simulation from the initial state.

So maybe not committing the solvers to one unique initialization recipe and letting the users develop their own mix is fine? 🧙 🤣

Yes, the initialization procedure will have to be more sophisticated in the presence of more complex physics (like faults and fractures for instance). If necessary, each solver could have an associated Task to perform some physics-specific initialization tasks. However, letting users develop their own mix may be challenging, because it requires some understanding of the internals of GEOSX, I think.

[TotoGaz]

Maybe the grain cannot be the solver: you may have different combination of solvers, values, etc. that may require some different initializations.

[joshua-white]

I think the workflow is a good one, though I would emphasize that Step 4 may involve some nuance. If any stress or traction conditions lie outside the yield surface after Step 3, we would have to move the yield surface at these points to ensure everything is in static equilibrium. Otherwise there will be a sudden deformation at t=0.

The alternative is to do a geologic timescale simulation with inelasticity enabled, to gradually compute an in situ stress and inelastic state that makes the most sense. I think both options are available using the currrent implementation. The second will work already, though Step 3 above would now take several nonlinear timesteps. For the first, we just have to add a yield surface reset algorithm to the solid models (maybe 15 lines of code).

[klevzoff]

For Step 3, have you considered instead forming and solving a statically condensed system with a regular single-physics linear solver? I.e.

delta x_uu = -J_uu^-1 (r_u - J_up x_p - J_us x_s)

This would require writing a special variation of multiphase poromechanics FEM kernel that assembles the modified right-hand side. But to me it still seems like a more straightforward and efficient approach, although I might be missing something?

[francoishamon]

Yes, you are right, currently what I have implemented in PR #1927 is a simple implementation that uses MGR to do the static condensation, but if we decide to use this initialization approach, we may want to do something better along the lines of what you suggest.

[joshua-white]

@francoishamon, was it not sufficient to just set dt = a big number in the standard assembly routine without any modifications? I would think this would automatically create the equilibrium problem we want to solve.

@klevzoff, I could see what you propose being more broadly useful in the context of one-way poromechanics (solve fluid flow on its own, and then postprocess for deformations).

[francoishamon]

@joshua-white
For our large-scale cases, I am facing two practical issues with the method that uses a coupled step with dt = a big number:

• For the properties that were given to us, this coupled initialization step is very difficult to converge for dt larger than 10 years because the nonlinear solver struggles. But maybe this time step size is enough to initialize the model, and maybe I could also improve the behavior of the solver by providing a better initial guess for the stress (better than 0 should not be too difficult)?
• The other practical issue is that I have to impose the datum pressure during the coupled initialization step (using our machinery for Dirichlet BC) to make sure that the pressure does not drift too far away from the data given to us. Currently in my tests I have not done that, but that should not be too difficult either, just a little bit tedious.

(and maybe the reason why I see the pressure drifting during the initialization step is because my initialization time step is not large enough, so the two problems are still too strongly coupled?)

The approach that I have used is similar to what @klevzoff is proposing but just more hacky (using MGR to do the static condensation). This initialization is doing one-way poromechanics by computing the stresses from the input pressure. The downside is that the system may not at a true equilibrium when we start the simulation, but the advantages seem to be that:

• The initialization is really easy to set up and everything is done internally (no need to add BCs to preserve the datum pressure)
• It is fast and only requires an elastic solve.

I think there is value in having both methods, one being rigorous, the other being fast and easy to use. But I don't have enough experience to make a decision between the two approaches (for instance, how far from equilibrium is the second approach?)