bug fix for min/max values not working on discontinuous problem - #1226
bug fix for min/max values not working on discontinuous problem#1226ee-nn wants to merge 11 commits into
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1226 +/- ##
==========================================
+ Coverage 95.38% 95.64% +0.25%
==========================================
Files 53 53
Lines 4140 4202 +62
==========================================
+ Hits 3949 4019 +70
+ Misses 191 183 -8 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
| # the extrema exports read the solution on the submesh of the volume | ||
| # subdomain, so check that the species is actually defined there. | ||
| # a field that is not a Species with a list of subdomains (a bare name, | ||
| # or a species generated from a trap) is skipped: it fails earlier, for | ||
| # unrelated reasons | ||
| is_extremum = isinstance( | ||
| export, | ||
| exports.MaximumVolume | ||
| | exports.MinimumVolume | ||
| | exports.MaximumSurface | ||
| | exports.MinimumSurface, | ||
| ) | ||
| if is_extremum and isinstance(export.field, _species.Species): | ||
| if isinstance(export, exports.SurfaceQuantity): | ||
| volume = self.surface_to_volume[export.surface] | ||
| location = f"surface {export.surface.id}" | ||
| else: | ||
| volume = export.volume | ||
| location = f"volume {volume.id}" | ||
| subdomains = export.field.subdomains | ||
| if isinstance(subdomains, list) and volume not in subdomains: | ||
| raise ValueError( | ||
| f"Cannot compute {export.title}: species " | ||
| f"{export.field.name} is not defined in the volume subdomain " | ||
| f"{volume.id} that {location} belongs to" | ||
| ) | ||
|
|
There was a problem hiding this comment.
This just provides some error catching with nice-looking raises so could be deleted if desired
There was a problem hiding this comment.
I think it's good to catch this early
| bot_min = F.MinimumVolume( | ||
| field=H, volume=bottom_volume, filename=f"{tmpdir}/bot_min.csv" | ||
| ) | ||
| bot_max = F.MaximumVolume(field=H, volume=bottom_volume) | ||
| top_min = F.MinimumVolume(field=H, volume=top_volume) | ||
| top_max = F.MaximumVolume(field=H, volume=top_volume) | ||
| bot_surf_min = F.MinimumSurface( | ||
| field=H, surface=bottom_surface, filename=f"{tmpdir}/bot_surf_min.csv" | ||
| ) | ||
| bot_surf_max = F.MaximumSurface(field=H, surface=bottom_surface) | ||
| top_surf_min = F.MinimumSurface(field=H, surface=top_surface) | ||
| top_surf_max = F.MaximumSurface(field=H, surface=top_surface) |
There was a problem hiding this comment.
Might be a little redundant, perhaps only the bottom domain could be checked?
There was a problem hiding this comment.
I think it's ok to look out for everything
RemDelaporteMathurin
left a comment
There was a problem hiding this comment.
Thanks for this @ee-nn !!! A first pass of comments before I go more in depth
| meshtags = self.facet_meshtags if facet_meshtags is None else facet_meshtags | ||
|
|
||
| if meshtags is None: | ||
| raise ValueError( |
There was a problem hiding this comment.
Could be an assertion instead with a message if false
| values = solution.x.array[dofs] | ||
|
|
||
| # a process may hold no dof of the surface at all, np.max would then raise | ||
| local_max = np.max(values) if values.size > 0 else -np.inf |
| # the extrema exports read the solution on the submesh of the volume | ||
| # subdomain, so check that the species is actually defined there. | ||
| # a field that is not a Species with a list of subdomains (a bare name, | ||
| # or a species generated from a trap) is skipped: it fails earlier, for | ||
| # unrelated reasons | ||
| is_extremum = isinstance( | ||
| export, | ||
| exports.MaximumVolume | ||
| | exports.MinimumVolume | ||
| | exports.MaximumSurface | ||
| | exports.MinimumSurface, | ||
| ) | ||
| if is_extremum and isinstance(export.field, _species.Species): | ||
| if isinstance(export, exports.SurfaceQuantity): | ||
| volume = self.surface_to_volume[export.surface] | ||
| location = f"surface {export.surface.id}" | ||
| else: | ||
| volume = export.volume | ||
| location = f"volume {volume.id}" | ||
| subdomains = export.field.subdomains | ||
| if isinstance(subdomains, list) and volume not in subdomains: | ||
| raise ValueError( | ||
| f"Cannot compute {export.title}: species " | ||
| f"{export.field.name} is not defined in the volume subdomain " | ||
| f"{volume.id} that {location} belongs to" | ||
| ) | ||
|
|
There was a problem hiding this comment.
I think it's good to catch this early
| bot_min = F.MinimumVolume( | ||
| field=H, volume=bottom_volume, filename=f"{tmpdir}/bot_min.csv" | ||
| ) | ||
| bot_max = F.MaximumVolume(field=H, volume=bottom_volume) | ||
| top_min = F.MinimumVolume(field=H, volume=top_volume) | ||
| top_max = F.MaximumVolume(field=H, volume=top_volume) | ||
| bot_surf_min = F.MinimumSurface( | ||
| field=H, surface=bottom_surface, filename=f"{tmpdir}/bot_surf_min.csv" | ||
| ) | ||
| bot_surf_max = F.MaximumSurface(field=H, surface=bottom_surface) | ||
| top_surf_min = F.MinimumSurface(field=H, surface=top_surface) | ||
| top_surf_max = F.MaximumSurface(field=H, surface=top_surface) |
There was a problem hiding this comment.
I think it's ok to look out for everything
| ) | ||
|
|
||
| with pytest.raises(ValueError, match="facet meshtags are required"): | ||
| my_export.compute() |
There was a problem hiding this comment.
If meshtags are always required why are they not a mandatory argument of compute?
There was a problem hiding this comment.
Good point, I can think about this more and update
There was a problem hiding this comment.
Pushed an implementation, let me know what you think
This reverts commit 6f92175.
| u: dolfinx.fem.Function | None = None, | ||
| volume_meshtags: dolfinx.mesh.MeshTags | None = None, |
There was a problem hiding this comment.
For a given problem, one would never invoke compute() with different volume mesthags and different us right? I think they should be attributes instead and set inside problem.initialise() somewhere. (I think we have a initialise_exports method)
| defined on. This is the expected behaviour for the discontinuous | ||
| problem, where each volume subdomain owns its own submesh. | ||
| """ | ||
| solution = self.field.post_processing_solution if u is None else u |
There was a problem hiding this comment.
this could be a property
There was a problem hiding this comment.
or a method rather
| problem, where each volume subdomain owns its own submesh. | ||
| """ | ||
| solution = self.field.post_processing_solution if u is None else u | ||
| meshtags = self.volume_meshtags if volume_meshtags is None else volume_meshtags |
There was a problem hiding this comment.
this should be an attribute
There was a problem hiding this comment.
Does this only apply to the volume_meshtag class and then we keep the facet_meshtags argument for the surface?
There was a problem hiding this comment.
i think this would apply to both Volume quantities and Surface Quantities (if that's what you're asking)
There was a problem hiding this comment.
Ah yeah that's what I was getting at, thanks 🙂
| defined on. This is the expected behaviour for the discontinuous | ||
| problem, where each volume subdomain owns its own submesh. | ||
| """ | ||
| solution = self.field.post_processing_solution if u is None else u |
There was a problem hiding this comment.
at the beginning of compute we should have an assertion to make sure that meshtags have been given
There was a problem hiding this comment.
This would only be needed for when we are not on a submesh right? If we have a property like self.is_submesh as you suggest for the if meshtags is None statement, then this would go in the else branch of the conditional block, no?
There was a problem hiding this comment.
I guess so yes!
| ) | ||
|
|
||
| self.value = mesh.comm.allreduce(np.min(solution.x.array[dofs]), op=MPI.MIN) | ||
| if meshtags is None: |
There was a problem hiding this comment.
this should maybe be a property too like self.is_submesh.
Again, one day we'll get rid of the distinction between discontinous and non-discontinuous #1191
Description
Provides a fix for issue #1225.
Summary
Essentially #1225 was the result of two separate problems:
initialise_exports()ofHydrogenTransportProblemDiscontinuousoverrides the base method fromHydrogenTransportProblem, but it never initializes its ownexport.volume_meshtagsorexport.facet_meshtags. Since it doesn't invokesuper()either, these attributes ofexportaren't carried over to the discontinuous class.Even with the meshtags attached,
compute()in Min/MaxVolume and Min/MaxSurface readself.field.post_processing_solution, which staysNonein the discontinuous case because the solution lives per-submesh inspecies.subdomain_to_post_processing_solution[subdomain].Motivation and Context
Type of Change
Testing
pytest)Specifically, new system test for the discontinuous case (which is essentially the MWE listed in #1225 - not sure if this is overkill?), and some smaller assertion tests to ensure that
meshtagsare always defined.Code Quality Checklist
ruff format .)ruff check .)Documentation