diff --git a/demos/cahnhilliard/demo_cahnhilliard.py.rst b/demos/cahnhilliard/demo_cahnhilliard.py.rst deleted file mode 100644 index 658efcf..0000000 --- a/demos/cahnhilliard/demo_cahnhilliard.py.rst +++ /dev/null @@ -1,158 +0,0 @@ -Cahn-Hilliard Demo -================== - -We are reprising the Cahn-Hilliard demo using :math:`C^1` Bell elements from -Kirby, Mitchell, "Code generation for generally mapped finite -elements," ACM TOMS 45(4), 2019. But now, we use Irksome to perform -the time stepping. - -The model is, with :math:`\Omega` the unit square, - -.. math:: - - \frac{\partial c}{\partial t} - \nabla \cdot M \left(\nabla\left(\frac{d f}{d c} - - \lambda \nabla^{2}c\right)\right) = 0 \quad {\rm in} - \ \Omega. - -where :math:`f` is some typically non-convex function -(in our case, we take :math:`f(c) = 100c^2(1-c)^2`, and -:math:`\lambda` and :math:`M` are scalar parameters controlling -rates. We are considering only the simple case of :math:`M` constant, -independent of :math:`c`. - -We close the system with boundary conditions - -.. math:: - - M\left(\nabla\left(\frac{d f}{d c} - \lambda \nabla^{2}c\right)\right) - \cdot n &= 0 \quad {\rm on} \ \Gamma, - - M \lambda \nabla c \cdot n &= 0 \quad {\rm on} \ \Gamma - -For simplicity, we'll use a direct solver at each time step. - -Boilerplate imports:: - - from firedrake import * - from firedrake.pyplot import tripcolor - import numpy as np - import os - from irksome import Dt, GaussLegendre, MeshConstant, TimeStepper - -We create a directory to store some output pictures:: - - if not os.path.exists("pictures"): - os.makedirs("pictures") - elif not os.path.isdir("pictures"): - raise RuntimeError("Cannot create output directory, file of given name exists") - -Set up the mesh and approximating space, including some refined ones -to allow visualizing our higher-order element on a :math:`P^1` space:: - - N = 16 - msh = UnitSquareMesh(N, N) - - vizmesh = MeshHierarchy(msh, 2)[-1] - - V = FunctionSpace(msh, "Bell", 5) - -Cahn-Hilliard parameters:: - - lmbda = Constant(1.e-2) - M = Constant(1) - - - def dfdc(cc): - return 200*(cc*(1-cc)**2-cc**2*(1-cc)) - -With Bell elements, the path of least resistance for strong boundary -conditions is a Nitsche-type method. Here is the parameter:: - - beta = Constant(250.0) - -Set up the time variables and a seeded initial condition:: - - MC = MeshConstant(msh) - dt = MC.Constant(5.0e-6) - T = MC.Constant(5.0e-6) - t = MC.Constant(0.0) - - np.random.seed(42) - c = Function(V) - c.dat.data[::6] = 0.63 + 0.2*(0.5 - np.random.random(c.dat.data[::6].shape)) - - v = TestFunction(V) - -Now we define the semidiscrete variational problem:: - - def lap(u): - return div(grad(u)) - - n = FacetNormal(msh) - h = CellSize(msh) - - F = (inner(Dt(c), v) * dx + - inner(M*grad(dfdc(c)), grad(v))*dx + - inner(M*lmbda*lap(c), lap(v))*dx - - inner(M*lmbda*lap(c), dot(grad(v), n))*ds - - inner(M*lmbda*dot(grad(c), n), lap(v))*ds + - inner(beta/h*M*lmbda*dot(grad(c), n), dot(grad(v), n))*ds) - -Bell elements are fourth-order accurate in :math:`L^2`, so we'll use a -time-stepping scheme of comparable accuracy:: - - butcher_tableau = GaussLegendre(2) - -Because of the nonlinear problem, we'll need to set set some Newton -parameters as well as the linear solver:: - - params = {'snes_monitor': None, 'snes_max_it': 100, - 'snes_linesearch_type': 'l2', - 'ksp_type': 'preonly', - 'pc_type': 'lu', 'mat_type': 'aij', - 'pc_factor_mat_solver_type': 'mumps'} - -Set up the output:: - - output = Function(FunctionSpace(vizmesh, "P", 1), - name="concentration") - - P5 = Function(FunctionSpace(msh, "P", 5)) - intp = Interpolator(c, P5) - - def interpolate_output(): - intp.interpolate() - return prolong(P5, output) - -Save the initial condition to a file:: - - import matplotlib.pyplot as plt - interpolate_output() - cs = tripcolor(output, vmin=0, vmax=1) - plt.colorbar(cs) - plt.savefig('pictures/init.pdf', format='pdf', bbox_inches='tight', pad_inches=0) - -Now let's set up the time stepper:: - - stepper = TimeStepper(F, butcher_tableau, t, dt, c, - solver_parameters=params) - -And advance the solution in time:: - - while float(t) < float(T): - if (float(t) + float(dt)) >= 1.0: - dt.assign(1.0 - float(t)) - stepper.advance() - t.assign(float(t) + float(dt)) - print(float(t), float(dt)) - -We'll save a snapshout of the final state:: - - interpolate_output() - cs = tripcolor(output, vmin=0, vmax=1) - plt.colorbar(cs) - plt.savefig('pictures/final.pdf', format='pdf', bbox_inches='tight', pad_inches=0) - -And report the amount of overshoot we get in the method:: - - print(np.max(c.dat.data[::6])) diff --git a/demos/monodomain/demo_monodomain_FHN.py.rst b/demos/monodomain/demo_monodomain_FHN.py.rst index daadc31..1625451 100644 --- a/demos/monodomain/demo_monodomain_FHN.py.rst +++ b/demos/monodomain/demo_monodomain_FHN.py.rst @@ -34,7 +34,7 @@ We start with standard Firedrake/Irksome imports:: And we set up the mesh and function space. Note this demo uses serendipity elements, but could just as easily use Lagrange on quads or triangles.:: - mesh = RectangleMesh(100, 100, 70, 70, quadrilateral=True) + mesh = RectangleMesh(20, 20, 70, 70, quadrilateral=True) polyOrder = 2 V = FunctionSpace(mesh, "S", polyOrder) diff --git a/docs/source/index.rst b/docs/source/index.rst index 5d167b9..541f017 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -64,13 +64,6 @@ to deploy more efficient methods: demos/demo_heat_rana.py demos/demo_dirk_parameters.py -Yes, Irksome can handle nonlinear problems as well: - -.. toctree:: - :maxdepth: 1 - - demos/demo_cahnhilliard.py - We now have support for DIRKs: .. toctree:: diff --git a/tests/test_split.py b/tests/test_split.py index 679a7d7..cd8973c 100644 --- a/tests/test_split.py +++ b/tests/test_split.py @@ -7,6 +7,7 @@ errornorm, grad, inner, pi, sin, split) from irksome import Dt, MeshConstant, RadauIIA, TimeStepper from irksome.tools import AI, IA +from pyop2.mpi import COMM_WORLD @pytest.mark.parametrize('splitting', (AI, IA)) @@ -124,7 +125,7 @@ def NavierStokesSplitTest(N, num_stages, Fimp, Fexp): MC = MeshConstant(mesh) t = MC.Constant(0.0) - dt = MC.Constant(0.5 / N) + dt = MC.Constant(1.0 / N) z_imp = Function(Z) z_split = Function(Z) @@ -138,22 +139,26 @@ def Ffull(z, test): bcs = [DirichletBC(Z.sub(0), as_vector([x*(1-x), 0]), (4,)), DirichletBC(Z.sub(0), as_vector([0, 0]), (1, 2, 3))] - nsp = [(1, VectorSpaceBasis(constant=True))] + nsp = [(1, VectorSpaceBasis(constant=True, comm=COMM_WORLD))] lunl = { - "mat_type": "aij", + "snes_type": "newtonls", "snes_rtol": 1e-10, "snes_atol": 1e-10, + "snes_converged_reason": None, + "snes_linesearch_type": "l2", "ksp_type": "preonly", "pc_type": "lu", - "pc_factor_shift_type": "nonzero"} + "pc_factor_mat_solver_type": "mumps" + } lulin = { "mat_type": "aij", "snes_type": "ksponly", "ksp_type": "preonly", "pc_type": "lu", - "pc_factor_shift_type": "nonzero"} + "pc_factor_mat_solver_type": "pastix", + "pc_factor_shift_type": "inblocks"} F_full = Ffull(z_imp, test_z) F_imp = Fimp(z_split, test_z) @@ -174,9 +179,8 @@ def Ffull(z, test): num_its_initial=10, num_its_per_step=4) - while (float(t) < 1.0): - if (float(t) + float(dt) > 1.0): - dt.assign(1.0 - float(t)) + nsteps = 6 + for _ in range(nsteps): imp_stepper.advance() imex_stepper.advance() t.assign(float(t) + float(dt)) @@ -196,7 +200,7 @@ def Ffull(z, test): def test_SplitNavierStokes(N, num_stages, Fimp, Fexp): error = NavierStokesSplitTest(N, num_stages, Fimp, Fexp) print(abs(error)) - assert abs(error) < 4e-8 + assert abs(error) < 4e-7 if __name__ == "__main__": diff --git a/tests/test_stokes.py b/tests/test_stokes.py index 452a03e..d721664 100644 --- a/tests/test_stokes.py +++ b/tests/test_stokes.py @@ -11,7 +11,11 @@ def StokesTest(N, butcher_tableau, stage_type="deriv", splitting=AI): - mesh = UnitSquareMesh(N, N) + mesh0 = UnitSquareMesh(N//4, N//4) + mh = MeshHierarchy(mesh0, 2) + mesh = mh[-1] + + ns = butcher_tableau.num_stages Ve = VectorElement("CG", mesh.ufl_cell(), 2) Pe = FiniteElement("CG", mesh.ufl_cell(), 1) @@ -45,22 +49,42 @@ def StokesTest(N, butcher_tableau, stage_type="deriv", splitting=AI): u, p = z.subfunctions u.interpolate(uexact) - - lu = {"mat_type": "aij", - "snes_type": "newtonls", - "snes_linesearch_type": "l2", - "snes_linesearch_monitor": None, - "snes_monitor": None, - "snes_rtol": 1e-8, - "snes_atol": 1e-8, - "snes_force_iteration": 1, - "ksp_type": "preonly", - "pc_type": "lu", - "pc_factor_mat_solver_type": "mumps"} + p.interpolate(pexact-assemble(pexact*dx)) + + ns = butcher_tableau.num_stages + ind_pressure = ",".join([str(2*i+1) for i in range(ns)]) + solver_params = { + "mat_type": "aij", + "snes_type": "ksponly", + "ksp_type": "fgmres", + "ksp_max_it": 200, + "ksp_gmres_restart": 30, + "ksp_rtol": 1.e-8, + "ksp_atol": 1.e-13, + "pc_type": "mg", + "pc_mg_type": "multiplicative", + "pc_mg_cycles": "v", + "mg_levels": { + "ksp_type": "chebyshev", + "ksp_max_it": 3, + "ksp_convergence_test": "skip", + "pc_type": "python", + "pc_python_type": "firedrake.ASMVankaPC", + "pc_vanka_construct_dim": 0, + "pc_vanka_sub_sub_pc_type": "lu", + "pc_vanka_sub_sub_pc_factor_mat_solver_type": "umfpack", + "pc_vanka_exclude_subspaces": ind_pressure}, + "mg_coarse": { + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + "mat_mumps_icntl_14": 200} + } stepper = TimeStepper(F, butcher_tableau, t, dt, z, stage_type=stage_type, - bcs=bcs, solver_parameters=lu, nullspace=nsp) + bcs=bcs, solver_parameters=solver_params, + nullspace=nsp) while (float(t) < 1.0): if (float(t) + float(dt) > 1.0): @@ -81,8 +105,10 @@ def StokesTest(N, butcher_tableau, stage_type="deriv", splitting=AI): # and check that the velocity is the right size (as observed from # a "by-hand" backward Euler code in Firedrake def NSETest(butch, stage_type, splitting): - N = 16 + N = 4 M = UnitSquareMesh(N, N) + mh = MeshHierarchy(M, 2) + M = mh[-1] V = VectorFunctionSpace(M, "CG", 2) W = FunctionSpace(M, "CG", 1) @@ -106,15 +132,39 @@ def NSETest(butch, stage_type, splitting): nullspace = [(1, VectorSpaceBasis(constant=True))] - solver_parameters = {"mat_type": "aij", - "snes_type": "newtonls", - "snes_linesearch_type": "bt", - "snes_rtol": 1e-10, - "snes_atol": 1e-10, - "snes_force_iteration": 1, - "ksp_type": "preonly", - "pc_type": "lu", - "pc_factor_mat_solver_type": "mumps"} + ns = butch.num_stages + ind_pressure = ",".join([str(2*i+1) for i in range(ns)]) + solver_params = { + "mat_type": "aij", + "snes_type": "newtonls", + "snes_converged_reason": None, + "snes_linesearch_type": "l2", + "snes_rtol": 1.e-8, + "ksp_rtol": 1.e-10, + "ksp_atol": 1.e-13, + "ksp_type": "fgmres", + "ksp_monitor": None, + "ksp_max_it": 200, + "ksp_gmres_restart": 30, + "pc_type": "mg", + "pc_mg_type": "multiplicative", + "pc_mg_cycles": "v", + "mg_levels": { + "ksp_type": "chebyshev", + "ksp_max_it": 3, + "ksp_convergence_test": "skip", + "pc_type": "python", + "pc_python_type": "firedrake.ASMVankaPC", + "pc_vanka_construct_dim": 0, + "pc_vanka_sub_sub_pc_type": "lu", + "pc_vanka_sub_sub_pc_factor_mat_solver_type": "umfpack", + "pc_vanka_exclude_subspaces": ind_pressure}, + "mg_coarse": { + "ksp_type": "preonly", + "pc_type": "lu", + "pc_factor_mat_solver_type": "mumps", + "mat_mumps_icntl_14": 200} + } MC = MeshConstant(M) t = MC.Constant(0.0) @@ -123,7 +173,7 @@ def NSETest(butch, stage_type, splitting): t, dt, up, bcs=bcs, stage_type="value", - solver_parameters=solver_parameters, + solver_parameters=solver_params, nullspace=nullspace) tfinal = 1.0 @@ -141,12 +191,12 @@ def NSETest(butch, stage_type, splitting): @pytest.mark.parametrize('stage_type', ("deriv", "value")) @pytest.mark.parametrize('splitting', (AI, IA)) -@pytest.mark.parametrize('N', [2**j for j in range(2, 4)]) +@pytest.mark.parametrize('N', [2**j for j in range(3, 4)]) @pytest.mark.parametrize('time_stages', (2, 3)) @pytest.mark.parametrize('butch', (LobattoIIIC, RadauIIA)) def test_Stokes(N, butch, time_stages, stage_type, splitting): error = StokesTest(N, butch(time_stages), stage_type, splitting) - assert abs(error) < 4e-9 + assert abs(error) < 1e-8 @pytest.mark.parametrize('stage_type', ("deriv", "value")) @@ -158,4 +208,5 @@ def test_NSE(butch, time_stages, stage_type): if __name__ == "__main__": - test_Stokes(8, RadauIIA, 2, "value", IA) + # test_Stokes(4, RadauIIA, 2, "deriv", IA) + test_NSE(RadauIIA, 2, "deriv")