GEOS
GEOS copied to clipboard
GEOS-DEV
Fixed issue #2853
This PR fixed issue #2853 .
The nonlinear system solver (solveNonlinearSystem) may invoke the updateState function, which can modify the stencil weights at each newton iteration. In this PR, the stencil weights are updated only once in implicitStepComplete function at the end of each solver step.
The issue is that this will effectively change the solution. I can see how it would help with convergence because it removes a nonlinearity but it is essentially treating the main source of coupling explicitly.
The issue is that this will effectively change the solution. I can see how it would help with convergence because it removes a nonlinearity but it is essentially treating the main source of coupling explicitly.
To fully account for the nonlinearity, we need to include the derivatives of m_weights w.r.t displacement and pressure, since it is a function of aperture. However, I think this is not essential.
https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L483-L484
The updating scheme of m_weights in SurfaceElementStencil.hpp is lagging, just as stated by @karimifard in https://github.com/GEOS-DEV/GEOS/issues/2853#issuecomment-1856491359
Essentially the update of the aperture in the flow discretization is lagging and the Jacobian accounts only for the impact of the aperture on permeability variation.
Yes, it is lagging but it's happening within the nonlinear iterations so when the solution has converged at a given time step all quantities are up to date.
The issue is that this will effectively change the solution. I can see how it would help with convergence because it removes a nonlinearity but it is essentially treating the main source of coupling explicitly.
To fully account for the nonlinearity, we need to include the derivatives of
m_weightsw.r.t displacement and pressure, since it is a function ofaperture. However, I think this is not essential. https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L483C1-L484C89
We do account for derivatives of the transmissibility w.r.t. displacement and pressure in the flux terms. The only thing that was done was to lag part of the update by 1 iteration but, as pointed out by Mohammad, at the end of the solution procedure for a timestep everything is up to date.
The updating scheme of
w_weightsinSurfaceElementStencil.hppis lagging, just as stated by @karimifard in #2853 (comment)Essentially the update of the aperture in the flow discretization is lagging and the Jacobian accounts only for the impact of the aperture on permeability variation.
Updating the stencil weights at the end of each Newton iteration is a standard Newton iteration, not lagging. $$J(y_k^{n+1})(y_{k+1}^{n+1} - y_k^{n+1}) = -R(y_k^{n+1})$$
$k$ is the Newton iteration count and $n$ is the time increment count.
If we try to use the standard Newton iteration, we should account for the contribution of the derivatives of m_weights (not stencil weights) in dt0_dvar1, dt1_dvar1, dt0_dvar2, and dt1_dvar2 as well, since the m_weights is a function of aperture.
https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L483-L484
https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L497-L498
Updating the stencil weights at the end of each Newton iteration is a standard Newton iteration, not lagging. J(ykn+1)(yk+1n+1−ykn+1)=−R(ykn+1)
k is the Newton iteration count and n is the time increment count.
If we try to use the standard Newton iteration, we should account for the contribution of the derivatives of
m_weights(not stencilweights) indt0_dvar1,dt1_dvar1,dt0_dvar2, anddt1_dvar2as well, since them_weightsis a function ofaperture.https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L483-L484
https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L497-L498
Yes, we are all well aware of this. The choice was made for a few reasons.
My initial point, though, was that this PR goes in the complete opposite direction by introducing an explicit update of the dependency.
Updating the stencil weights at the end of each Newton iteration is a standard Newton iteration, not lagging. J(ykn+1)(yk+1n+1−ykn+1)=−R(ykn+1) k is the Newton iteration count and n is the time increment count. If we try to use the standard Newton iteration, we should account for the contribution of the derivatives of
m_weights(not stencilweights) indt0_dvar1,dt1_dvar1,dt0_dvar2, anddt1_dvar2as well, since them_weightsis a function ofaperture. https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L483-L484https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L497-L498
Yes, we are all well aware of this. The choice was made for a few reasons.
My initial point, though, was that this PR goes in the complete opposite direction by introducing an explicit update of the dependency.
Yes. In this PR, the surface element stencil weights is updated by old hydraulic aperture aperture0, while the transmissibility component $\frac{w^2}{12}$ is still updated by current hydraulic aperture aperture.
Updating the stencil weights at the end of each Newton iteration is a standard Newton iteration, not lagging. J(ykn+1)(yk+1n+1−ykn+1)=−R(ykn+1) k is the Newton iteration count and n is the time increment count. If we try to use the standard Newton iteration, we should account for the contribution of the derivatives of
m_weights(not stencilweights) indt0_dvar1,dt1_dvar1,dt0_dvar2, anddt1_dvar2as well, since them_weightsis a function ofaperture. https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L483-L484https://github.com/GEOS-DEV/GEOS/blob/decc28110acf58053a0d6911f5ea7a4d9d7f76d1/src/coreComponents/finiteVolume/SurfaceElementStencil.hpp#L497-L498
Yes, we are all well aware of this. The choice was made for a few reasons. My initial point, though, was that this PR goes in the complete opposite direction by introducing an explicit update of the dependency.
Yes. In this PR, the surface element stencil
weightsis updated by old hydraulic apertureaperture0, while the transmissibility component w212 is still updated by current hydraulic apertureaperture.
If it really solves convergence issues for some of those test cases we can maybe keep it as optional by introducing a flag.
I have tested the pknViscosityDominated_benchmark.xml with this PR on my side and the results match well.
The MPI parameters I used for the case is as follows:
export OMP_NUM_THREADS=2
mpirun -np 45 geos -i ./pknViscosityDominated_benchmark.xml -x 3 -y 3 -z 5
Following figure is the results from this PR hydrofractureSolver and it took about 01h13m04s to complete the run
Following figure is the results from hydrofractureSolver before merging PR #2623 and it took about 01h44m03s to complete the run
hydrofractureSolver from the latest develop branch has convergence issues (performance degradation #2853 ) and requires frequent time cutback for most time increment.
Please note that I changed the penaltyStiffness value from 1.0e0 to 1.0e12, otherwise the simulation would stop around 60s.
https://github.com/GEOS-DEV/GEOS/blob/741b6ab459c01266627b518df67769dd0dbf6700/inputFiles/hydraulicFracturing/pknViscosityDominated_base.xml#L55
Any updates on this PR?
No further updates. I believe this Pull Request should resolve the performance problem. If required, I can modify this Pull Request to address any merge conflicts.
Any updates on this PR?
No further updates. I believe this Pull Request should resolve the performance problem. If required, I can modify this Pull Request to address any merge conflicts.
Since it is going to be hard to agree what way is the right, how about we make an option (flag) to switch the way the weights are updated. @CusiniM @karimifard what do you think?
Any updates on this PR?
No further updates. I believe this Pull Request should resolve the performance problem. If required, I can modify this Pull Request to address any merge conflicts.
Since it is going to be hard to agree what way is the right, how about we make an option (flag) to switch the way the weights are updated. @CusiniM @karimifard what do you think?
I am fine if we have a flag to choose how the stencils are updated.
@liyangrock are you willing to implement a flag to switch between current and new way ?
@liyangrock are you willing to implement a flag to switch between current and new way ?
Hi @paveltomin , I originally thought that you were planning to adopt new weight coefficients, so I didn't update this PR. I will implement a flag to switch between the two methods.
Hi @paveltomin, I’ve submitted a new PR to tackle issue #2853, please see #3108.