Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Global treatment of out of balance forces. #2633

Merged
merged 5 commits into from Oct 25, 2019
Merged

Conversation

@endJunction
Copy link
Member

endJunction commented Aug 27, 2019

Follows #2635, where the <initial_stress> tag is added. Reenables the tests commented in #2635.

  1. Feature description was added to the changelog
  2. Tests covering your feature were added?
  3. Any new feature or behavior change was documented?
@endJunction endJunction force-pushed the endJunction:FoobGlobal branch 2 times, most recently from cf29062 to a9693cf Aug 27, 2019
@@ -184,6 +184,17 @@ void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
// equation every time and could not forget it.
}

void NonlinearSolver<NonlinearSolverTag::Newton>::calculateOutOfBalanceForces(

This comment has been minimized.

Copy link
@norihiro-w

norihiro-w Aug 27, 2019

Collaborator

IMO, NumLib should be independent of physics.

This comment has been minimized.

Copy link
@endJunction

endJunction Aug 28, 2019

Author Member

I just updated the documentation, there are some more details.
In the document we named the initial residuum F_oob -- out-of-balance forces. These are not the "physical" forces, but general external loads. For mechanics, these are indeed forces as

$F^oob = r_0 = \int\Omega B^T \sigma_0\; d\Omega$

This comment has been minimized.

Copy link
@norihiro-w

norihiro-w Aug 29, 2019

Collaborator

Thanks for the doc. Is there any way to see the documentation without doing git lfs?

Do we need this calculation besides mechanics? To me, "load", "force" sounds specific to mechanics.

This comment has been minimized.

Copy link
@endJunction

endJunction Aug 29, 2019

Author Member

the doc.

https://deploy-preview-2633--ogs.netlify.com/docs/benchmarks/small-deformations/mechanics-linear-nonequilibrium-states/

The preview to the benchmark sites is under the "deploy/netlify -- Details" in the checks of the PR. There in "docs", "Selected Benchmarks", you will find the above site.

"load", "force" ...

Same arguments are valid for pressure or saturation. One could rename Foob to initial_residuum, r_0, or something similar to avoid mechanics relation...

This comment has been minimized.

Copy link
@nagelt

nagelt Aug 29, 2019

Member

@endJunction The F^oob equation above also contains body forces and Neumann bc's right? I think you were just brief but want to make sure...

@norihiro-w I sent you the equations a while ago to replace the non-equ. stresses with a better more general solution. Indeed "forces" points to mechanics, like much of the FE nomenclature (e.g. stiffness matrix). So I wouldn't have a problem with that name. The most general alternative could be "out of balance fluxes". That applies to all processes. In case of mechanics, forces are linear momentum fluxes.... then we could even keep the F ;-)

This comment has been minimized.

Copy link
@nagelt

nagelt Aug 29, 2019

Member

and yes, it is intended that this applies to all processes where needed ... we can write this for a general balance equation and Foob would always be the initial out of balance flux that a user may want to drive an initial process or not (and hence would have to compensate) ...
time for dinner now in China ...

This comment has been minimized.

Copy link
@norihiro-w

norihiro-w Aug 29, 2019

Collaborator

@endJunction thank you for the tip. I could read it.

@endJunction @nagelt Indeed I'm not suggesting renaming the function. I'm suggesting moving this functionality to somewhere in ProcessLib. NumLib should be a general pure numerical library and we should avoid contaminate it with physics except a concept of space and time. I guess you can use a post-assembly function caller in the nonlinear solver class without changing the general interface. If you still want to change NumLib, it would be better if you use a concept in math/numeric, such as residual shift/offset (i'm not sure this example is appropriate).

@nagelt Is "Foob" a word you defined? I was asking Dima about the name because I'm surprised that this trivial algorithm has a name.

This comment has been minimized.

Copy link
@nagelt

nagelt Aug 30, 2019

Member

@norihiro-w maybe I misunderstand you and say something totally obvious to you now, but we just use "Foob" as a short-hand for $F_oob$ = out-of-balance-forces/fluxes. We could also call this initial residual if we leave physics out of the name ... but for most engineers the physical name implies the concept better than the purely mathematical one..

regarding ProcessLib: if we can move it there without much extra code into a general place, I agree. I'm not fully aware of the structure there and the effort involved, so I'll leave this call to @endJunction

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

I agree with @norihiro-w, this function should be in ProcessLib for the mechanical related processes.

@norihiro-w

This comment has been minimized.

Copy link
Collaborator

norihiro-w commented Aug 27, 2019

what is "the global Foob's procedure"? If there is any paper about it, would you please give me citation?

Copy link
Member

TomFischer left a comment

Looks technical good.

@endJunction endJunction force-pushed the endJunction:FoobGlobal branch from 9c7209d to 65feb6e Aug 28, 2019
@endJunction

This comment has been minimized.

Copy link
Member Author

endJunction commented Aug 28, 2019

what is "the global Foob's procedure"?

"global" is just in the sense that the implementation is valid for all nonlinear problems, where the residuum needs to be adjusted, whereas (there is an open PR #2561 being superseded by #2635 and another one) non-equilibrium initial state treatment happened on the ("local") process scale.

@endJunction endJunction mentioned this pull request Aug 28, 2019
3 of 3 tasks complete
@endJunction endJunction force-pushed the endJunction:FoobGlobal branch 2 times, most recently from 7739374 to 1d0394d Aug 29, 2019
@norihiro-w

This comment has been minimized.

Copy link
Collaborator

norihiro-w commented Aug 29, 2019

@nagelt why do you have to remove "r0" also from the mass balance equation in HM/RM? There should be no problem to calculate hydrostatic pressure. Maybe you have some applications demanding this, but I don't think it is usual case. It would be nicer if you introduce an option to activate or deactivate calculating r0 for non-mechanic process.

@nagelt

This comment has been minimized.

Copy link
Member

nagelt commented Aug 30, 2019

@norihiro-w It was a feature request to have this also for the mass balance. You're right, we want separate control over the individual process residuals, so that you can decide to only compensate for individual processes, e.g. mechanics. @endJunction will do this in a future PR.

@@ -184,6 +184,17 @@ void NonlinearSolver<NonlinearSolverTag::Newton>::assemble(
// equation every time and could not forget it.
}

void NonlinearSolver<NonlinearSolverTag::Newton>::calculateOutOfBalanceForces(

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

I agree with @norihiro-w, this function should be in ProcessLib for the mechanical related processes.

@@ -232,6 +243,10 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
sys.getJacobian(J);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());

// Apply out-of-balance forces if set.

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

It can be moved to assembleWithJacobianConcreteProcess.

This comment has been minimized.

Copy link
@endJunction

endJunction Sep 1, 2019

Author Member

Probably I'm missing something here, but I think it is not possible to move this part out of the nonlinear solver. The reason is the sys.getResidual(x, res) few lines above. Computation of the residuum depends on the time discretization scheme. Passing the time discretization scheme to the processes is not a good option.
The residuum can not be computed in the process in general case.

@norihiro-w @wenqing I'd like to see some sketches/ideas how to move the computation into the processes.

This comment has been minimized.

Copy link
@norihiro-w

norihiro-w Sep 3, 2019

Collaborator

@endJunction At least, for mechanics, you don't need to care the time discretization, because it's solving static problems. For dynamic problems, I have no idea how you can make F^oob consistent with any time discretization scheme.

return;

_equation_system->assemble(x);
_f_oob = &NumLib::GlobalVectorProvider::provider.getVector();

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

You can move _f_oob to class TimeLoop, and initialize it on demanding.

@@ -108,11 +110,18 @@ class NonlinearSolver<NonlinearSolverTag::Newton> final

void assemble(GlobalVector const& x) const override;

void calculateOutOfBalanceForces(GlobalVector const& x) override;

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

It is better to move it to Process. In future, we can use CRTP idiom to encapsulate all mechanical process related stuff.

@@ -127,11 +136,18 @@ class NonlinearSolver<NonlinearSolverTag::Newton> final
//! conservative approach.
double const _damping;

GlobalVector* _f_oob = nullptr; //!< Out-of-balance forces vector.

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

Better to move it to Process.

/// before the first time step. The forces are zero if the external forces
/// are in equilibrium with the initial state/initial conditions. During the
/// simulation the new residuum reads \f$ \tilde r = r - F_{\rm oob} \f$.
bool _compensate_initial_out_of_balance_forces = false;

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

It it better to move it to Process.

NonlinearSolverStatus solve(
GlobalVector& x,
std::function<void(int, GlobalVector const&)> const&
postIterationCallback) override;

void compensateInitialOutOfBalanceForces(bool const value)

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

It it better to move it Process.

@@ -167,6 +167,34 @@ std::vector<GlobalVector*> setInitialConditions(
return process_solutions;
}

void calculateOutOfBalanceForces(

This comment has been minimized.

Copy link
@wenqing

wenqing Aug 30, 2019

Member

can be a member of Process

@endJunction endJunction force-pushed the endJunction:FoobGlobal branch 2 times, most recently from 7ebc220 to e500ebb Aug 30, 2019
@endJunction

This comment has been minimized.

Copy link
Member Author

endJunction commented Aug 30, 2019

Just rebased. I'll try to move the code to ProcessLib/Process. I'm not sure whether this is good idea, because as I understood it, the procedure is independent of any process implementations and is a nonlinear solver feature:

r_0 = F(u_0)

solve

F(u) - r_0 = 0;

where the original formulation is to minimize F(u) is changed to minimize F(u) - F(u_0).

Furthermore, the non-equilibrium conditions can be given through pressure/temperature or saturation, and not only for mechanical processes.

@endJunction endJunction force-pushed the endJunction:FoobGlobal branch 3 times, most recently from 898b283 to 2130f0c Aug 31, 2019
@wenqing
wenqing approved these changes Sep 5, 2019
Copy link
Member

wenqing left a comment

From the point view of numerics, the nonlinear solver should not included this part of computation. We can consider it again once the non-linear part is refactored.

@endJunction endJunction added the WIP 🏗 label Sep 5, 2019
@endJunction

This comment has been minimized.

Copy link
Member Author

endJunction commented Sep 5, 2019

From the point view of numerics, the nonlinear solver should not included this part of computation.

With this I need some help to understand it; the nonlinear solver is, until now, solving an equation of this type:

F(u) = 0

for u, using either fixed point iterations or some Newton scheme. Here we modify the equation to be solved to

F(u) - F(u_0) = 0

Now, I'm not against it if this can be moved to ProcessLib but I admit not to have an idea how to.

@ThieJan

This comment has been minimized.

Copy link
Contributor

ThieJan commented Sep 6, 2019

I do not have the complete overview, but if its neither process lib nor nonlinear solver, maybe we have to add something new in between. E.g., wrap around the class responsible for the evaluation of F(u) returning the 'shifted' value?

Copy link
Member

wenqing left a comment

After the discussion with @endJunction, I think that the current place for the additional initial RHS is adequate. The only suggestions are the naming of balance forces, and adding the documentation about the new nonlinear equation of F(u) -F(u_0)=0

@endJunction endJunction force-pushed the endJunction:FoobGlobal branch 2 times, most recently from 9847f16 to 451baf9 Oct 21, 2019
@endJunction

This comment has been minimized.

Copy link
Member Author

endJunction commented Oct 23, 2019

Rebased. Renamed the "out-of-balance forces" to "non-equilibrium initial residuum".

@endJunction endJunction force-pushed the endJunction:FoobGlobal branch 2 times, most recently from ce1271f to 1953ff3 Oct 23, 2019
@endJunction endJunction added please review and removed WIP 🏗 labels Oct 23, 2019
@endJunction endJunction force-pushed the endJunction:FoobGlobal branch from 1953ff3 to b442962 Oct 24, 2019
@TomFischer TomFischer merged commit 77676d2 into ufz:master Oct 25, 2019
3 checks passed
3 checks passed
continuous-integration/jenkins/pr-merge This commit looks good
Details
deploy/netlify Deploy preview ready!
Details
ufz.ogs #20191024.2 succeeded
Details
@TomFischer TomFischer deleted the endJunction:FoobGlobal branch Oct 25, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
6 participants
You can’t perform that action at this time.