From f611804208e8586935a1d96c6d18e74c7496369c Mon Sep 17 00:00:00 2001 From: gdmcbain Date: Tue, 17 Dec 2019 14:40:03 +1100 Subject: [PATCH] default AMG-preconditioned solvers from pyamgcl #5 --- .../st08_navier_stokes_cylinder.py | 21 ++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/08_navier_stokes_cylinder/st08_navier_stokes_cylinder.py b/08_navier_stokes_cylinder/st08_navier_stokes_cylinder.py index b9dfb92..fee3dd7 100644 --- a/08_navier_stokes_cylinder/st08_navier_stokes_cylinder.py +++ b/08_navier_stokes_cylinder/st08_navier_stokes_cylinder.py @@ -12,6 +12,8 @@ from skfem.models.general import divergence from skfem.models.poisson import laplace, vector_laplace +import pyamgcl + @skfem.bilinear_form def vector_mass(u, du, v, dv, w): @@ -99,6 +101,16 @@ def embed(xy: np.ndarray) -> np.ndarray: return np.pad(xy, ((0, 0), (0, 1)), 'constant') +solvers = [pyamgcl.solver(pyamgcl.amg(skfem.condense(K_lhs, + D=dirichlet['u'], + expand=False))), + pyamgcl.solver(pyamgcl.amg(skfem.condense(L['p'], + D=dirichlet['p'], + expand=False))), + pyamgcl.solver(pyamgcl.amg(skfem.condense(M / dt, + D=dirichlet['u'], + expand=False)))] + with TimeSeriesWriter(Path(__file__).with_suffix('.xdmf').name) as writer: writer.write_points_cells(embed(mesh.p.T), {'triangle': mesh.t.T}) @@ -112,17 +124,20 @@ def embed(xy: np.ndarray) -> np.ndarray: uv = skfem.solve(*skfem.condense( K_lhs, K_rhs @ uv_ - P @ (2 * p_ - p__) - skfem.asm(acceleration, basis['u'], w=basis['u'].interpolate(u)), - uv0, D=dirichlet['u'])) + uv0, D=dirichlet['u']), + solver=solvers[0]) # Step 2: pressure correction dp = skfem.solve(*skfem.condense(L['p'], (B / dt) @ uv, - D=dirichlet['p'])) + D=dirichlet['p']), + solver=solvers[1]) # Step 3: velocity correction p = p_ + dp - du = skfem.solve(*skfem.condense(M / dt, -P @ dp, D=dirichlet['u'])) + du = skfem.solve(*skfem.condense(M / dt, -P @ dp, D=dirichlet['u']), + solver=solvers[2]) u = uv + du uv_ = uv