Skip to content

Commit

Permalink
Merge pull request #43 from pyroteus/23_remove_class_attrs
Browse files Browse the repository at this point in the history
Avoid class attributes
  • Loading branch information
jwallwork23 authored Jan 15, 2024
2 parents 2062f12 + 56fde15 commit fee6497
Showing 1 changed file with 28 additions and 14 deletions.
42 changes: 28 additions & 14 deletions movement/monge_ampere.py
Original file line number Diff line number Diff line change
Expand Up @@ -48,9 +48,6 @@ class MongeAmpereMover_Base(PrimeMover):
Base class for mesh movers based on the solution
of Monge-Ampere type equations.
"""
residual_l2_form = 0
norm_l2_form = 0
phi_old = None

def __init__(self, mesh, monitor_function, **kwargs):
"""
Expand Down Expand Up @@ -106,6 +103,10 @@ def apply_initial_guess(self, phi_init=None, sigma_init=None, **kwargs):
:kwarg sigma_init: initial guess for the Hessian
"""
if phi_init is not None and sigma_init is not None:
assert hasattr(self, "phi")
assert hasattr(self, "sigma")
assert hasattr(self, "phi_old")
assert hasattr(self, "sigma_old")
self.phi.assign(phi_init)
self.sigma.assign(sigma_init)
self.phi_old.assign(phi_init)
Expand All @@ -114,10 +115,10 @@ def apply_initial_guess(self, phi_init=None, sigma_init=None, **kwargs):
raise ValueError("Need to initialise both phi *and* sigma")

@property
@PETSc.Log.EventDecorator("MongeAmpereBase.diagnostics")
def diagnostics(self):
@PETSc.Log.EventDecorator()
def _diagnostics(self):
"""
Compute diagnostics:
Compute the following diagnostics:
1) the ratio of the smallest and largest element volumes;
2) equidistribution of elemental volumes;
3) relative L2 norm residual.
Expand All @@ -129,8 +130,10 @@ def diagnostics(self):
w *= w
std = np.sqrt(w.sum()/w.size)
equi = std/mean
residual_l2 = firedrake.assemble(self.residual_l2_form).dat.norm
norm_l2 = firedrake.assemble(self.norm_l2_form).dat.norm
assert hasattr(self, "_residual_l2_form")
assert hasattr(self, "_norm_l2_form")
residual_l2 = firedrake.assemble(self._residual_l2_form).dat.norm
norm_l2 = firedrake.assemble(self._norm_l2_form).dat.norm
residual_l2_rel = residual_l2/norm_l2
return minmax, residual_l2_rel, equi

Expand Down Expand Up @@ -159,6 +162,7 @@ def l2_projector(self):
"""
if hasattr(self, '_l2_projector'):
return self._l2_projector
assert hasattr(self, "phi_old")
u_cts = firedrake.TrialFunction(self.P1_vec)
v_cts = firedrake.TestFunction(self.P1_vec)

Expand Down Expand Up @@ -266,8 +270,8 @@ def __init__(self, mesh, monitor_function, **kwargs):
self.theta_form = self.monitor*ufl.det(I + self.sigma_old)*self.dx
self.residual = self.monitor*ufl.det(I + self.sigma_old) - self.theta
psi = firedrake.TestFunction(self.P1)
self.residual_l2_form = psi*self.residual*self.dx
self.norm_l2_form = psi*self.theta*self.dx
self._residual_l2_form = psi*self.residual*self.dx
self._norm_l2_form = psi*self.theta*self.dx

@property
@PETSc.Log.EventDecorator("MongeAmpereMover.create_pseudotimestepper")
Expand All @@ -277,6 +281,8 @@ def pseudotimestepper(self):
"""
if hasattr(self, '_pseudotimestepper'):
return self._pseudotimestepper
assert hasattr(self, "phi")
assert hasattr(self, "phi_old")
phi = firedrake.TrialFunction(self.P1)
psi = firedrake.TestFunction(self.P1)
a = ufl.inner(ufl.grad(psi), ufl.grad(phi))*self.dx
Expand All @@ -302,6 +308,8 @@ def equidistributor(self):
"""
if hasattr(self, '_equidistributor'):
return self._equidistributor
assert hasattr(self, "phi")
assert hasattr(self, "sigma")
if self.dim != 2:
raise NotImplementedError # TODO
n = ufl.FacetNormal(self.mesh)
Expand All @@ -324,6 +332,10 @@ def move(self):
"""
Run the relaxation method to convergence and update the mesh.
"""
assert hasattr(self, "phi")
assert hasattr(self, "sigma")
assert hasattr(self, "phi_old")
assert hasattr(self, "sigma_old")
for i in range(self.maxiter):

# L2 project
Expand All @@ -342,7 +354,7 @@ def move(self):
self.theta.assign(firedrake.assemble(self.theta_form) / self.total_volume)

# Check convergence criteria
minmax, residual, equi = self.diagnostics
minmax, residual, equi = self._diagnostics
if i == 0:
initial_norm = residual
PETSc.Sys.Print(f"{i:4d}"
Expand Down Expand Up @@ -404,8 +416,8 @@ def __init__(self, mesh, monitor_function, **kwargs):
self.theta_form = self.monitor*ufl.det(I + self.sigma_old)*self.dx
self.residual = self.monitor*ufl.det(I + self.sigma_old) - self.theta
psi = firedrake.TestFunction(self.P1)
self.residual_l2_form = psi*self.residual*self.dx
self.norm_l2_form = psi*self.theta*self.dx
self._residual_l2_form = psi*self.residual*self.dx
self._norm_l2_form = psi*self.theta*self.dx

@property
@PETSc.Log.EventDecorator("MongeAmpereMover.create_equidistributor")
Expand All @@ -415,6 +427,7 @@ def equidistributor(self):
"""
if hasattr(self, '_equidistributor'):
return self._equidistributor
assert hasattr(self, "phisigma")
if self.dim != 2:
raise NotImplementedError # TODO
n = ufl.FacetNormal(self.mesh)
Expand All @@ -432,6 +445,7 @@ def update_monitor(cursol):
"""
Callback for updating the monitor function.
"""
assert hasattr(self, "phisigma_old")
with self.phisigma_old.dat.vec as v:
cursol.copy(v)
self.l2_projector.solve()
Expand Down Expand Up @@ -472,7 +486,7 @@ def monitor(snes, i, rnorm):
firedrake.assemble(self.L_P0, tensor=self.volume)
self.volume.interpolate(self.volume/self.original_volume)
self.mesh.coordinates.assign(self.xi)
minmax, residual, equi = self.diagnostics
minmax, residual, equi = self._diagnostics
PETSc.Sys.Print(f"{i:4d}"
f" Min/Max {minmax:10.4e}"
f" Residual {residual:10.4e}"
Expand Down

0 comments on commit fee6497

Please sign in to comment.