Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Average surface and volume derived quantities in Cylindrical and Spherical #820

Open
wants to merge 18 commits into
base: main
Choose a base branch
from

Conversation

jhdark
Copy link
Collaborator

@jhdark jhdark commented Jul 26, 2024

Proposed changes

Describe the big picture of your changes here to communicate to the maintainers why we should accept this pull request. If it fixes a bug or resolves a feature request, be sure to link to that issue.

Types of changes

What types of changes does your code introduce to FESTIM?

  • Bugfix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)
  • Breaking change (fix or feature that would cause existing functionality to not work as expected)
  • Code refactoring
  • Documentation Update (if none of the other choices apply)
  • New tests

Checklist

  • Black formatted
  • Unit tests pass locally with my changes
  • I have added tests that prove my fix is effective or that my feature works
  • I have added necessary documentation (if appropriate)

Further comments

If this is a relatively large or complex change, kick off the discussion by explaining why you chose the solution you did and what alternatives you considered, etc...

Copy link

codecov bot commented Jul 26, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 99.57%. Comparing base (282071e) to head (e747d52).
Report is 26 commits behind head on main.

Additional details and impacted files
@@           Coverage Diff           @@
##             main     #820   +/-   ##
=======================================
  Coverage   99.56%   99.57%           
=======================================
  Files          61       61           
  Lines        2747     2796   +49     
=======================================
+ Hits         2735     2784   +49     
  Misses         12       12           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@RemDelaporteMathurin RemDelaporteMathurin added the enhancement New feature or request label Jul 29, 2024
@jhdark jhdark marked this pull request as ready for review August 5, 2024 14:05
Copy link
Collaborator

@RemDelaporteMathurin RemDelaporteMathurin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here are a few comments.

On a side note, this makes me think we should try and test for cases where the surface varies in both $r$ and $z$. Like this:

image

@@ -47,3 +48,101 @@ def compute(self):
return f.assemble(self.function * self.ds(self.surface)) / f.assemble(
1 * self.ds(self.surface)
)


class AverageSurfaceCylindrical(AverageSurface):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since #832 you need to override the allowed_meshes property for these quantities.

Is this where we could make use of things like abstract methods and stuff?

Comment on lines +140 to +142
# dV_z = r dr dtheta , assuming axisymmetry dV_z = theta r dr
# dV_r = r dz dtheta , assuming axisymmetry dV_r = theta r dz
# in both cases the expression with self.dx is the same
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not true for spherical, please adapt

return avg_surf


class AverageSurfaceSpherical(AverageSurface):
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since spherical coordinates can only be used in 1D, then the surface is always 0D.

Do we need this class at all? Since for 0D surfaces $\int f \ dS = f$ then:

$$ \frac{\int f \ ds }{\int 1 \ ds} = \frac{f}{1} = f $$

Then we have AverageSurfaceSpherical = AverageSurface.

Does that make sense?

top_id = 2
top_surface.mark(surface_markers, top_id)

my_export = AverageSurfaceCylindrical("solute", top_id)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These tests pass if you replace the class on this line with AverageSurface (cartesian).

Maybe the test case isn't chosen correctly chosen

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's because the concentration field is constant on the top surface...

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Try with

c_fun = lambda r: 2 * r

right_surface.mark(surface_markers, right_id)
ds = f.Measure("ds", domain=mesh_fenics, subdomain_data=surface_markers)

my_export = AverageSurfaceSpherical("solute", right_id)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Follwing my comment above (AverageSurfaceSpherical = AverageSurface) if you replace this by AverageSurface the tests pass

volume_markers.set_all(1)
dx = f.Measure("dx", domain=mesh_fenics, subdomain_data=volume_markers)

my_export = AverageVolumeCylindrical("solute", 1)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment as above. If you change this to simply AverageVolume then the tests pass. The c_fun function doesn't vary in $r$ so you're not really testing that your implementaiton works well here

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Try with $c = 3r$


my_export = AverageVolumeSpherical("solute", 1)
V = f.FunctionSpace(mesh_fenics, "P", 1)
c_fun = lambda r: c_left + (c_right - c_left) / (r1 - r0) * r
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just noting that $c$ is not equal to $c_\mathrm{right}$ with this expression (replace $r$ by $r1$)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants