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.

Sign up for GitHub

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

Copy mesh upon creating Mover #48

Merged
merged 3 commits into from
Jan 25, 2024
Merged

Copy mesh upon creating Mover #48

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0314926
#45: Copy mesh upon creating Mover
jwallwork23 Jan 23, 2024
88aae71
#45: It doesn't make sense to create a new Mover on adapted mesh
jwallwork23 Jan 23, 2024
b04f422
#45: Ensure tests use mover.mesh post-adapt
jwallwork23 Jan 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions movement/monge_ampere.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,17 +77,17 @@ def __init__(self, mesh, monitor_function, **kwargs):
super().__init__(mesh, monitor_function=monitor_function)

# Create function spaces
self.P0 = firedrake.FunctionSpace(mesh, "DG", 0)
self.P1 = firedrake.FunctionSpace(mesh, "CG", 1)
self.P1_vec = firedrake.VectorFunctionSpace(mesh, "CG", 1)
self.P1_ten = firedrake.TensorFunctionSpace(mesh, "CG", 1)
self.P0 = firedrake.FunctionSpace(self.mesh, "DG", 0)
self.P1 = firedrake.FunctionSpace(self.mesh, "CG", 1)
self.P1_vec = firedrake.VectorFunctionSpace(self.mesh, "CG", 1)
self.P1_ten = firedrake.TensorFunctionSpace(self.mesh, "CG", 1)

# Create objects used during the mesh movement
self.theta = firedrake.Constant(0.0)
self.monitor = firedrake.Function(self.P1, name="Monitor function")
self.monitor.interpolate(self.monitor_function(self.mesh))
self.volume = firedrake.Function(self.P0, name="Mesh volume")
self.volume.interpolate(ufl.CellVolume(mesh))
self.volume.interpolate(ufl.CellVolume(self.mesh))
self.original_volume = firedrake.Function(self.volume)
self.total_volume = firedrake.assemble(firedrake.Constant(1.0) * self.dx)
self.L_P0 = firedrake.TestFunction(self.P0) * self.monitor * self.dx
Expand Down
20 changes: 11 additions & 9 deletions movement/mover.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,28 +12,30 @@ class PrimeMover(object):
"""

def __init__(self, mesh, monitor_function=None, **kwargs):
self.mesh = mesh
self.mesh = firedrake.Mesh(mesh.coordinates.copy(deepcopy=True))
self.monitor_function = monitor_function
self.dim = mesh.topological_dimension()
self.dim = self.mesh.topological_dimension()
if self.dim != 2:
raise NotImplementedError(
f"Dimension {self.dim} has not been considered yet"
)
self.gdim = mesh.geometric_dimension()
self.gdim = self.mesh.geometric_dimension()
self.plex = self.mesh.topology_dm
self.vertex_indices = self.plex.getDepthStratum(0)
self.edge_indices = self.plex.getDepthStratum(1)

# Measures
degree = kwargs.get("quadrature_degree")
self.dx = firedrake.dx(domain=mesh, degree=degree)
self.ds = firedrake.ds(domain=mesh, degree=degree)
self.dS = firedrake.dS(domain=mesh, degree=degree)
self.dx = firedrake.dx(domain=self.mesh, degree=degree)
self.ds = firedrake.ds(domain=self.mesh, degree=degree)
self.dS = firedrake.dS(domain=self.mesh, degree=degree)

# Mesh coordinate functions
self.coord_space = mesh.coordinates.function_space()
self._x = firedrake.Function(mesh.coordinates, name="Physical coordinates")
self.xi = firedrake.Function(mesh.coordinates, name="Computational coordinates")
self.coord_space = self.mesh.coordinates.function_space()
self._x = firedrake.Function(self.mesh.coordinates, name="Physical coordinates")
self.xi = firedrake.Function(
self.mesh.coordinates, name="Computational coordinates"
)
self.v = firedrake.Function(self.coord_space, name="Mesh velocity")

def _get_coordinate_section(self):
Expand Down
6 changes: 3 additions & 3 deletions movement/spring.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,10 +39,10 @@ def __init__(self, mesh, **kwargs):
:arg mesh: the physical mesh
"""
super().__init__(mesh)
self.HDivTrace = firedrake.FunctionSpace(mesh, "HDiv Trace", 0)
self.HDivTrace_vec = firedrake.VectorFunctionSpace(mesh, "HDiv Trace", 0)
self.HDivTrace = firedrake.FunctionSpace(self.mesh, "HDiv Trace", 0)
self.HDivTrace_vec = firedrake.VectorFunctionSpace(self.mesh, "HDiv Trace", 0)
self.f = firedrake.Function(self.mesh.coordinates.function_space())
self.displacement = np.zeros(mesh.num_vertices())
self.displacement = np.zeros(self.mesh.num_vertices())

@property
@PETSc.Log.EventDecorator("SpringMover_Base.facet_areas")
Expand Down
19 changes: 8 additions & 11 deletions test/test_monge_ampere.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ def test_uniform_monitor(method, exports=False):
mover = MongeAmpereMover(mesh, const_monitor, method=method)
num_iterations = mover.move()

assert np.allclose(coords, mesh.coordinates.dat.data)
assert np.allclose(coords, mover.mesh.coordinates.dat.data)
assert num_iterations == 0


Expand All @@ -42,14 +42,11 @@ def test_continue(method, exports=False):
# Solve the problem to a weak tolerance
mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=0.1)
num_it_init = mover.move()
phi, sigma = mover.phi, mover.sigma
if exports:
File("outputs/init.pvd").write(phi, sigma)
File("outputs/init.pvd").write(mover.phi, mover.sigma)

# Continue with a new Mover
mover = MongeAmpereMover(
mesh, ring_monitor, method=method, rtol=rtol, phi_init=phi, sigma_init=sigma
)
Comment on lines -49 to -52
Copy link
Member Author

Choose a reason for hiding this comment

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

It doesn't really make sense to 'continue' the mesh movement by creating a new Mover with the computational mesh being the adapted mesh from an earlier Mover. They have different baselines.

Copy link
Collaborator

Choose a reason for hiding this comment

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

Ah yes, I was wondering about that/discussing with Matt: in our current implementation we are always working with a uniform mesh as the baseline right? And so the monitor we provide is specifying the required variation in the output physical mesh regardless of the previous state of the mesh (when calling the mover multiple times). We could also imagine a different approach where the monitor expresses the variation with respect to the previous mesh, where if the mesh is already optimal the monitor would be constant.

# Continue with a tighter tolerance
mover.rtol = rtol
num_it_continue = mover.move()
if exports:
File("outputs/continue.pvd").write(mover.phi, mover.sigma)
Expand Down Expand Up @@ -83,14 +80,14 @@ def test_change_monitor(method, exports=False):
mover.move()
if exports:
File("outputs/ring.pvd").write(mover.phi, mover.sigma)
assert not np.allclose(coords, mesh.coordinates.dat.data, atol=tol)
assert not np.allclose(coords, mover.mesh.coordinates.dat.data, atol=tol)

# Adapt to a constant monitor
mover.monitor_function = const_monitor
mover.move()
if exports:
File("outputs/const.pvd").write(mover.phi, mover.sigma)
assert np.allclose(coords, mesh.coordinates.dat.data, atol=tol)
assert np.allclose(coords, mover.mesh.coordinates.dat.data, atol=tol)


@pytest.mark.slow
Expand All @@ -113,10 +110,10 @@ def test_bcs(method, fix_boundary):
mover.move()

# Check boundary lengths are preserved
bnd_new = assemble(one * ds(domain=mesh))
bnd_new = assemble(one * ds(domain=mover.mesh))
assert np.isclose(bnd, bnd_new)

# Check boundaries are indeed fixed
if fix_boundary:
bnd_coords_new = mesh.coordinates.dat.data[bnodes]
bnd_coords_new = mover.mesh.coordinates.dat.data[bnodes]
assert np.allclose(bnd_coords, bnd_coords_new)
Loading