Add function to PETSc interface to return initial condition problem #1443
Conversation
dallan-keylogic
requested review from
andrewlee94,
bpaul4 and
lbianchi-lbl
as code owners
|
This looks really useful, thanks for making the change. I'll give it a review soon. |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1443 +/- ##
==========================================
- Coverage 76.36% 76.36% -0.01%
==========================================
Files 394 394
Lines 64953 64993 +40
Branches 14404 14412 +8
==========================================
+ Hits 49601 49630 +29
- Misses 12793 12803 +10
- Partials 2559 2560 +1 ☔ View full report in Codecov by Sentry. |
|
@Robbybp, just a friendly nag to review when (or if) you have some time. |
ksbeattie
requested review from
adam-a-a
and removed request for
dangunter,
blnicho,
lbianchi-lbl and
bpaul4
|
@Robbybp , could you please review this sometime soon? I'd really like it to make the August release. |
| def get_initial_condition_problem( | ||
| model, | ||
| time, | ||
| initial_time, | ||
| representative_time=None, | ||
| initial_constraints=None, | ||
| initial_variables=None, | ||
| detect_initial=True, | ||
| flattened_problem=None, | ||
| ): |
There was a problem hiding this comment.
Is there any reason initial_time is required? I would think that time.first() would be a reasonable default.
| if flattened_problem is None: | ||
| if representative_time is None: | ||
| raise RuntimeError( | ||
| "The user must supply either the flattened problem or a representative time." | ||
| ) | ||
| flattened_problem = _get_flattened_problem( | ||
| model=model, time=time, representative_time=representative_time | ||
| ) |
There was a problem hiding this comment.
Why is either a flattened problem or representative time required? Shouldn't this function just use the same defaults that petsc_dae_by_time_element uses?
There was a problem hiding this comment.
I'm trying to avoid having to flatten the problem multiple times. I haven't profiled it, but I suspect it's rather slow. I debated reimplementing petsc_dae_by_time_element as an Object so that the flattened problem can be cached. Maybe that will happen eventually.
There was a problem hiding this comment.
I have no problem with a flattened_problem argument, I just think if neither it nor representative time is provided, we should default to _get_flattened_problem(model=model, time=time, representative_time=time.at(2)).
| detect_initial (bool): If True, add non-time-indexed variables and | ||
| constraints to initial_variables and initial_constraints. |
There was a problem hiding this comment.
I don't love the name of this option, but since it is the same as in petsc_dae_by_time_element, I think it should stay.
| flattened_problem (dict): Dictionary returned by get_flattened_problem. | ||
| If not provided, get_flattened_problem will be called at representative_time. |
There was a problem hiding this comment.
I assume this option is provided so we don't have the overhead of recomputing the flattened model when calling within petsc_dae_by_time_element? I'm fine with this option, but would prefer for it to not be part of this function's public API. This is mostly because I don't want a user to have to construct a dict with undocumented keys (or use the private _get_flattened_problem function) to use this option. I recommend renaming the option to _flattened_problem, and adding to the docstring that its behavior is subject to change.
For future applications where a "flattened model" is necessary, the DynamicModelInterface from pyomo.contrib.mpc may be useful (see https://github.com/Pyomo/pyomo/blob/404fd6d997d24f7209e8af4c427a13d6e10f19c4/pyomo/contrib/mpc/interfaces/model_interface.py#L52 or https://pyomo.readthedocs.io/en/stable/contributed_packages/mpc/interface.html#pyomo.contrib.mpc.interfaces.model_interface.DynamicModelInterface).
| if len(constraints) <= 0: | ||
| raise RuntimeError( | ||
| "Zero constraints in initial condition problem, therefore " | ||
| "there is no block to return." | ||
| ) |
There was a problem hiding this comment.
We could just return an empty block, or None, in this case. I think I like this better than raising an error here that we then catch above. We already check the number of constraints before trying to solve the block, so an empty block should be fine, although None might be more explicit.
| try: | ||
| init_subsystem = get_initial_condition_problem( | ||
| model=m, | ||
| time=time, | ||
| initial_time=t0, | ||
| initial_constraints=initial_constraints, | ||
| initial_variables=initial_variables, | ||
| detect_initial=detect_initial, | ||
| flattened_problem=flattened_problem, | ||
| ) | ||
| except RuntimeError as err: | ||
| if "Zero constraints" in err.message: | ||
| init_subsystem = None | ||
| else: | ||
| raise err |
There was a problem hiding this comment.
See below. I'd rather this function return None, and we handle it here, than rely on parsing the error message.
|
Progress not expected until end of Sept (per @dallan-keylogic) |
Fixes
NoneSummary/Motivation:
When using the PETSc interface, a square initial condition problem is first solved (with either PETSc-snes or IPOPT) to make sure the initial condition is consistent before sending integrations problems to PETSc-TS. However, when you end up with degree of freedom problems or a degenerate problem, it's difficult to debug. As an expedient, I would edit the PETSc problem to return the problem in the middle of the context manager.
I've added a function to get this problem directly and modified the PETSc interface to use this function to get the initial condition problem to solve to ensure consistency.
Legal Acknowledgement
By contributing to this software project, I agree to the following terms and conditions for my contribution: