Tupek/system solver extended#1609
Conversation
… we go along. workout some stress output details.
…eanup cycle zero handling.
| weak_form_->AddDomainIntegral( | ||
| Dimension<spatial_dim>{}, DependsOn<all_params...>{}, | ||
| [dt, cycle, mode, integrand](double time, auto X, auto... inputs) { | ||
| return integrand(TimeInfo(time, *dt, *cycle, *mode), X, inputs...); |
There was a problem hiding this comment.
I think reconstructing time_info here will also result in double adding dt, as we observed in time_discretized_weak_form. Since we are already storing dt_, cycle_, and mode_, can we store time_info_ = time_info when we invoke residual() and jacobian()?
Then the lambda injected for these add integral functions would be
[time_info, integrand] (double /* time */, auto X, auto... inputs) {
return integrand(time_info, X, inputs...);
}
| return applyTimeRuleToPrefix( | ||
| rule, t_info, | ||
| [&](auto... self_states_and_tail) { | ||
| constexpr std::size_t n_self = sizeof...(self_states_and_tail) - tail_count; | ||
| auto all = std::forward_as_tuple(self_states_and_tail...); | ||
| return [&]<std::size_t... Si, std::size_t... Ti>(std::index_sequence<Si...>, std::index_sequence<Ti...>) { | ||
| return applyCouplingTimeRules( | ||
| coupling, t_info, | ||
| [&](auto... interpolated_coupling) { | ||
| return std::forward<Callback>(callback)(std::get<Si>(all)..., interpolated_coupling...); | ||
| }, | ||
| std::get<n_self + Ti>(all)...); | ||
| }(std::make_index_sequence<n_self>{}, std::make_index_sequence<tail_count>{}); | ||
| }, | ||
| raw_args...); | ||
| } |
There was a problem hiding this comment.
The nested return here is not very human-readable. Can we simplify this method?
| smith::TimeInfo(physics->time(), dt, static_cast<size_t>(physics->cycle()))); | ||
| } | ||
|
|
||
| std::cout << "reaction norm: " << physics->getReactionStates().front().get()->Norml2() << '\n'; |
There was a problem hiding this comment.
I'd prefer SLIC_INFO + std::format, easier to read imo.
| int main(int argc, char* argv[]) | ||
| { | ||
| // _init_start | ||
| smith::ApplicationManager application_manager(argc, argv); |
There was a problem hiding this comment.
These examples are user facing so should have some inline comments
| #include "smith/differentiable_numerics/nonlinear_solve.hpp" | ||
| #include "smith/differentiable_numerics/make_time_info_material.hpp" | ||
| #include "smith/physics/functional_objective.hpp" | ||
| #include "gretl/wang_checkpoint_strategy.hpp" |
There was a problem hiding this comment.
I think this is a larger problem than this PR but we have gretl as an optional TPL guarded by the CMake option SMITH_ENABLE_GRETL. This will fail if someone turns that off, which I don't think anyone has tried.
There was a problem hiding this comment.
That's a good point. Gretl really isn't optional, core functionality (differentiability) depends on it. Unless I'm mistaken? If not, the CMake option should be removed IMO.
There was a problem hiding this comment.
Issue created. we can handle it outside of this PR #1630
btalamini
left a comment
There was a problem hiding this comment.
Creating multi-physics is very nice with this!
|
|
||
| /// @brief evaluate value of the ode state as used by the integration rule | ||
| template <typename T1> | ||
| /// @brief Return the static field value. |
There was a problem hiding this comment.
Looks like the agent added a second set of doxygen comments, (one before the template, the other after. Delete one set. The new set seems more specific then the old, maybe delete the old ones.
| /// @c h = 1e-4 * max(1,|t|). | ||
| /// | ||
| /// Rebuilt on each call from current value-level essentials, so late additions are reflected. | ||
| const smith::BoundaryConditionManager& getSecondDerivativeManager() const; |
There was a problem hiding this comment.
What is the purpose of this? It doesn't seem to be used. Is this to enable acceleration BCs?
| Strumpack, /**< Strumpack MPI-enabled direct frontal solver*/ | ||
| PetscCG, /**< PETSc MPI-enabled conjugate gradient solver */ | ||
| PetscGMRES, /**< PETSc MPI-enabled generalize minimal residual solver */ | ||
| None /**< Preconditioner application only, No linear solver Krylov iterations */ |
There was a problem hiding this comment.
Maybe the intent would be clearer to call this one PrecondOnly or something like that?
chapman39
left a comment
There was a problem hiding this comment.
It would be good to track this PRs new feature(s) in the changelog's unreleased section under "added"
https://github.com/llnl/smith/blob/develop/RELEASE-NOTES.md?plain=1#L22
| addBodyIntegralImpl(body_name, integrand, std::integer_sequence<int, active_parameters...>{}); | ||
| } | ||
|
|
||
| /// @brief Add a body integral to the weak form. |
There was a problem hiding this comment.
| /// @brief Add a body integral to the weak form. | |
| /// @brief Add a body integral depending on all input fields to the weak form. |
| addBoundaryIntegralImpl(boundary_name, integrand, std::integer_sequence<int, active_parameters...>{}); | ||
| } | ||
|
|
||
| /// @brief Add a boundary integral to the weak form. |
There was a problem hiding this comment.
| /// @brief Add a boundary integral to the weak form. | |
| /// @brief Add a boundary integral depending on all input fields to the weak form. |
| addInteriorBoundaryIntegralImpl(interior_name, integrand, std::integer_sequence<int, active_parameters...>{}); | ||
| } | ||
|
|
||
| /// @brief Add an interior boundary integral to the weak form. |
There was a problem hiding this comment.
| /// @brief Add an interior boundary integral to the weak form. | |
| /// @brief Add an interior boundary integral depending on all input fields to the weak form. |
| /// @brief Accessor to get a reference to the underlying ShapeAwareFunctional in case more direct access is needed. | ||
| /// @return Reference to ShapeAwareFunctional instance. | ||
| ShapeAwareFunctional<ShapeDispSpace, OutputSpace(InputSpaces...)>& getShapeAwareResidual() { return *weak_form_; } | ||
| protected: |
There was a problem hiding this comment.
Are these methods still necessary? Overloaded functions to add integrals with all inputs are defined previously already.
| addBodyIntegralImpl(body_name, qfunction, std::integer_sequence<int, active_parameters...>{}); | ||
| } | ||
|
|
||
| /// @brief Add a body integral to the objective function. |
There was a problem hiding this comment.
| /// @brief Add a body integral to the objective function. | |
| /// @brief Add a body integral depending on all input fields. |
| addBoundaryIntegralImpl(boundary_name, qfunction, std::integer_sequence<int, active_parameters...>{}); | ||
| } | ||
|
|
||
| /// @brief Add a boundary integral to the objective function. |
There was a problem hiding this comment.
| /// @brief Add a boundary integral to the objective function. | |
| /// @brief Add a boundary integral depending on all input fields. |
| { | ||
| dt_ = time_info.dt(); | ||
| cycle_ = time_info.cycle(); | ||
| current_time_info_ = &time_info; |
There was a problem hiding this comment.
Why do we have a guard to prevent current_time_info_ from being nullified prematurally in functional_weak_form but not for functional_objective
lihanyu97
left a comment
There was a problem hiding this comment.
The PR looks great and allows flexibility of constructing multiphysics modules from previously defined single physics modules, which can significantly speedup future developments.
Considering the substaintial changes, can you summarize a list of interface breaking changes (eg. removal of time_discretized_weak_form)? It will be easier for other PRs created before this one to adjust to the new API.
This introduces a new interface for constructing multiphysics systems of equations, their boundary conditions, time integration rules, and solving them using potentially customized staggered or coupled algorithms.