Reductor for the stationary Stokes equation#2448
Conversation
sdrave
left a comment
There was a problem hiding this comment.
Hi @ullmannsven, thanks for your work! I have only looked at the changes in pymor.algorithms so far. I think it would be good to put them into a separate PR to get parts of this merged faster. Also, we would have better documentation about these changes.
| zeros = [[None for _ in range(ncols)] for _ in range(nrows)] | ||
| for l in range(nrows): | ||
| for j in range(ncols): | ||
| if op.blocks[l, j] is not None: | ||
| source_ = op.blocks[l, j].source | ||
| range_ = op.blocks[l, j].range | ||
| zeros[l][j] = ZeroOperator(range_, source_) |
There was a problem hiding this comment.
You can pass None for a zero block to BlockOperator, so I don't think you need this.
There was a problem hiding this comment.
Actually, you do, because BlockOperator.__init__ will fail if have a row or column of Nones. However, it might make sense to give that method optional range and source space arguments to handle that special case. Then you could just pass op.range and op.source everywhere without having to worry about the zero operators.
There was a problem hiding this comment.
I followed this approach, it first led it to sum unexpected errors when a projected system got projected to a subbasis again. I found a fix, but i would like to here your opinion if this is now good enough.
Codecov Report❌ Patch coverage is
Additional details and impacted files
Flags with carried forward coverage won't be shown. Click here to find out more.
🚀 New features to boost your workflow:
|
|
@ullmannsven, could you rebase on main to get a clean commit history? |
d0cb99a to
a0ea1f2
Compare
|
I added a I therefore thought about introducing handling for blocked system for the |
I think the real issue here is that there is no clear concept of what should be in the I see two viable solutions:
I lean towards the second option. What do you think? |
I went with the second option. Is this what you had in mind? I’m basically only performing the projection when the first two assertions in project pass. |
sdrave
left a comment
There was a problem hiding this comment.
Hi @ullmannsven, sorry for the long delay!
| if dims == self._last_rom_dims: | ||
| return self._last_rom | ||
| else: | ||
| elif dims != self._last_rom_dims: | ||
| return self._reduce_to_subbasis(dims) | ||
|
|
||
| return self._last_rom | ||
|
|
There was a problem hiding this comment.
This is required in case supremizer enrichment is used. If reduce gets called the first time (without the dims argument), the the lengths of the current RBs are stored in the dims dict. Due to supremizer enrichment, the lengths of the RBs change. These new lengths are stored in self.last_rom_dims, as this gets computed after _reduce, where the actually projection of the operators and the supremizer enrichment takes place. In the old version, the check dims != self._last_rom_dims would always return True in case supremizer enrichment is used, which is certainly not desired as then projection to the subbasis would be performed. Does that make sense?
There was a problem hiding this comment.
Ah, now I understand the issue. However, the new version is also problematic:
- Assume we are using
StationaryRBReductor. I initialise it with / extend its basis to a basis of dimension N+1. Then I callreductor.reduce(N). With your code I would get a ROM of order N+1.
-ForStationaryRBStokesReductor,reduce_to_subasisprobably does a different thing than what the user would expect. For instance, if I initialize it with a u basis of dimension N and later call reduce with dimension N-1 for the u basis, I would get a ROM where one supremizer has been removed instead of a u basis vector.
I think the only proper solution here is that the supremizers need to be kept separate from the u basis. Then in project_operators/project_operators_to_subbasis one needs to assemble the actual projection basis by appending the right amount of supremizers and orthogonalize them against the other basis vectors.
Downsides:
- The resulting ROM will have a larger solution space dimension than the sums of the provided dimensions. (But this already happens with your code when I initialize with some bases and then call
reducewith the dimensions of these bases.) - One cannot use
project_to_subbasisto have a faster projection.
Upsides:
- The ROM will always be the supremizer-enrichment ROM corresponding to the provided RB dims.
- Each of the ROMs will be stable.
I don't think there are better alternatives, but maybe I'm mistaken?
|
|
||
|
|
||
| class StationaryLSRBReductor(ProjectionBasedReductor): | ||
| """Least-squares projection based reductor for stationary problems. |
There was a problem hiding this comment.
Is "least-squares projection" something people say? Does it really mean to minimize the norm of the resiual? Could this be misleading? (Maybe "least-squares Petrov-Galerkin projection" would be better?)
There was a problem hiding this comment.
Yes, maybe using Petrov-Galerkin projection here is a bit more clear (at least for the case without normal equations?) I dont have a better phrasing in mind at the moment, so i changed it to Least-squares Petrov-Galerkin. I also changed the next sentence in the documentation.
| 'output_functional': project_to_subbasis(rom.output_functional, None, dim_trial) | ||
| } | ||
|
|
||
| else: |
There was a problem hiding this comment.
I think you want to save the StationarayLSRBReductor as an attribute and call its project_operator_to_subasis here.
There was a problem hiding this comment.
I changed it, but it only works if project_operators_to_subbasis in the StationaryLSRBReductor gets an additional, optional argument last_rom (the StationaryLSRBReductor itself has no _last_rom, as reduce (where this is initialized) gets called on the StationaryRBStokesReductor by the user).
To be consistent with the other reductors in basic.py i changed this for all reductors. Is this okay? It should not break any other code and can be a nice option to have this available?
There was a problem hiding this comment.
Maybe it would then be better to call reduce on the sub-reductor and extract the operators from the resulting ROM?
Anyway, I just became aware the project_to_subbasis does not really do what it promises. By just adding the dimensions, you happen to truncate the last basis vectors, to whatever field they belong to. This is clearly not the intention of the user. I fear there is really no efficient implementation of project_operators_to_subbasis here. You just have to assemble a new effective basis an project again.
(Something could be done when the block structure would be preserved during projection. With supremizer enrichment, one would additionally have to do separate projections for the supremizer/non-supremizer parts of the basis and then glue everything together. I don't think it is worth the hassle to change everything just for a faster reduce_to_subbasis.)
|
Tests run fine now, but we get an
in the 'scikit_fem oldest' build. Scikit-FEM is pinned to 6.0.0 in the 'oldest' image. Could you check in which version |
|
This pull request modifies make ci_requirementsand commit the changed files to ensure that CI runs with the updated dependencies. |
…mplify.py and image.py
… the documentation
ff9010f to
214e15d
Compare
sdrave
left a comment
There was a problem hiding this comment.
This is, hopefully, my last round of review for this PR. Mostly minor comments. The only thing that might require relevant work is my comment regarding the discretization. But as I wrote: this can also wait until after the release.
New features
This PR introduces a
StationaryRBStokesReductorfor the stationary stokes equation.Implementation of three reduction methods: supremizer_galerkin, ls-ls (least squares and using least_square=True as solver option) and ls-normal (solves the least-squares problem via the normal equation).
Further adds the stationary stokes equation in 2D posed on a unit circle to model.examples.py, using skfem to discretize the problem and assemble the matrices (currently implemented via Taylor-Hood elements). One can choose the
rhsof the stokes equations, but not the domain or the discretization method at the moment.Open work: