idaholab
Adding direct acceleration calculation to central difference
This PR is addressing the first half of #28008, which to implement a method that allows for a solution update calculated by the direct calculation of acceleration from the internal and external forces.
Reason
This is the first step in improving explicit contact as described in #28008
Design
The "direct" solution update, is calculated as the following:
$a_n=M^{-1}(F^{ext}-F^{int})$ $v_{n+1/2}=v_{n-1/2}+\Delta t a_n$ $d_{n+1}=d_n + \Delta t v_{n+1/2}$
These steps can be found in ExplicitTimeIntegrator.C
The following changes have been made:
-
NonlinearSystemBase.h- Added
setNodalBCs()member function. This sets nodal BCs for use with the direct time integrator. It is identical tosetInitialSolution()but without predictor calls.
- Added
-
TimeIntegrator.h- Added function
isDirect()that returns a bool describing whether or not the direct time integrator is being used
- Added function
-
ExplicitTimeIntegrator.C- Added additional option for
is_directto compute the solution update from acceleration
- Added additional option for
-
ActuallyExplicitEuler.C- Added a call to
setNodalBCs()after solution update if using direct solver
- Added a call to
-
CentralDifference.C- Added a direct acceleration calculation if in the nonlinear system. This effectively does nothing as is, but will be useful for future work with explicit contact.
- Added parameter
use_directto allow a user to specify if they want to do use a direct acceleration calculation for time integration
-
InertialForce.C- Modified
computeQpResidual()to have an option to return 0 if using a direct time integrator - Modified
computeQpJacobian()to return the mass matrix without factor tags needed for the current method
- Modified
Impact
This will allow a user to specific use_direct=true to invoke a central difference time integrator that computes the solution update based on a direct calculation of acceleration from the internal and external forces.
@TheGreatCid are you working with anyone at the national labs on this? To know who to assign for the review
I'm working with Daniel Schwen. I did not realize drafts were public, this is not ready for review yet
Documentation still needs to be updated, however any code review would be appreciated
Job Documentation on e8072f8 wanted to post the following:
View the site here
This comment will be updated on new commits.
Job Coverage on e8072f8 wanted to post the following:
Framework coverage
| 5bf0d8 | #27991 e8072f | ||||
|---|---|---|---|---|---|
| Total | Total | +/- | New | ||
| Rate | 85.06% | 85.04% | -0.03% | 26.00% | |
| Hits | 104544 | 104552 | +8 | 13 | |
| Misses | 18360 | 18398 | +38 | 37 | |
Modules coverage
Solid mechanics
| 5bf0d8 | #27991 e8072f | ||||
|---|---|---|---|---|---|
| Total | Total | +/- | New | ||
| Rate | 84.83% | 84.83% | -0.00% | 75.00% | |
| Hits | 27780 | 27783 | +3 | 3 | |
| Misses | 4966 | 4967 | +1 | 1 | |
Full coverage reports
Reports
-
framework -
chemical_reactions -
combined -
contact -
electromagnetics -
external_petsc_solver -
fluid_properties -
fsi -
functional_expansion_tools -
geochemistry -
heat_transfer -
level_set -
misc -
navier_stokes -
optimization -
peridynamics -
phase_field -
porous_flow -
ray_tracing -
rdg -
reactor -
richards -
scalar_transport -
solid_mechanics -
solid_properties -
stochastic_tools -
thermal_hydraulics -
xfem
Warnings
-
frameworknew line coverage rate 26.00% is less than the suggested 90.0% -
solid_mechanicsnew line coverage rate 75.00% is less than the suggested 90.0%
This comment will be updated on new commits.
@lindsayad This implementation of CD really only works with solid mechanics due to how the mass matrix needs to be constructed. At least I am not familiar enough with fluid mechanics to know if the call to the mass matrix construction would return the correct matrix. Should I put in a warning to state this?
But I can also just ask: are you running into issues related to mass matrices and fluid mechanics? Or was this a preemptive question?
It's a preemptive question. My motivation for this was to fix explicit contact, so I did not give much though to other physics
My number one concern about this PR is whether it's too specific to solid mechanics such that all this work should be there instead of the framework. If this doesn't work for other physics, then maybe we are missing some generality in the framework that would enable you to implement the time integration you want in your module/application in a seamless way. I would ask you and @dschwen to weigh in and tell me whether you feel like you're locked into a box that doesn't fit very well and if that is the case, how we might make it better
We are certainly missing some generality here. You may be right that this would be best suited for the solid mechanics module. There may a way to make this more general, but I not familiar enough with the other modules to figure that out quickly.
The lack of generality can be seen in the changes I made to InertialForce.C. I need to modify the residual returned and ComputeQpJacobian as these had terms from the original CD formulation. I am assuming those exist in the other physics as well.
It might be overkill to change other time dependent kernels to might this formulation as the point of having this addition CD formulation is for explicit contact
The changes in InertialForce definitely look a bit funny as Jacobians are supposed to be the derivatives of the residual with respect to the solution dofs and you're returning a 0 residual ... yet the Jacobian is not 0. We introduced LinearSystem recently. I wonder if that would be a better paradigm for the work you're doing where you fill a RHS and a matrix (which would just be the mass matrix when using explicit time integration), and there is not the same J = dR/du conceptual relationship
Yea, I did know it was a bit of a bastardization of the method, but I wasn't aware of another option as that is how the mass matrix is returned in ActuallyExplicitEuler. I was not aware of LinearSystem, I'll need to take a look at that. It looks like I need to take another stab at this and restructure things. I appreciate the feedback!
Though the fact that the linear system does not support time integrators is an issue
Ha, yes that is a big issue. It would be great to have someone work on that, but it would make sense if that's not what you want to focus your attention on
If it didn't make sense or would take a lot of effort to re-use the current TimeIntegrator classes for LinearSystem, then this would be a great opportunity to use PETSc TS
Perhaps @grmnptr could weigh-in with his thoughts on time integration for LinearSystem
We talked about this in person and the conclusion was that it would be very hard to do this with the LinearSystem because at the moment it does not have support for:
- Finite Elements
- TimeIntegrators So a lot of work would be needed here. I recommended registering a new matrix to the nonlinear system in the time integrator and creating a custom ThreadedLoop to assemble the mass matrix into it. That might run into issues on the other fronts so feel free to share other suggestions!
If we are going that route, then I really think we should limit the number of changes made in the framework. I don't want to add second class support in the framework because we don't currently have the time to do the work that needs to be done for first class support. @recuero's explicit dynamics changes didn't really feel quite right when we merged them and I'm worried that this further jams the square peg into the round hole
The better route might be to move this time integrator into solid_mechanics, and then use the tagging system and MassMatrix.C to pass the mass matrix in. MassMatrix.C is something I was unaware of until yesterday evening.
@lindsayad I opened a new PR, https://github.com/idaholab/moose/pull/28175, with a different approach to implementing this time integrator
I won't review the other one yet as you have it marked as draft, but I do like it a lot better. It looks quite clean!
I left it as a draft because the documentation/tests are not yet implemented, but the code is ready for review. Should I mark it as ready now, or wait until I have added in doc/tests?
Might as well wait until you've added the documentation and tests so that a reviewer can look at it in one fell swoop. The input file can influence how I feel about the code design