Remove tests without taping from adjoint test directory. (#4870)

Author: JHopeCollins

tests/firedrake/adjoint/test_covariance_operator.py

@@ -1,235 +1,11 @@
 import pytest
-import numpy as np
-from scipy.sparse import csr_matrix
-import petsctools
 from firedrake import *
 from firedrake.adjoint import *
 
 
-def petsc2numpy_vec(petsc_vec):
-    """Allgather a PETSc.Vec."""
-    gvec = petsc_vec
-    gather, lvec = PETSc.Scatter().toAll(gvec)
-    gather(gvec, lvec, addv=PETSc.InsertMode.INSERT_VALUES)
-    return lvec.array_r.copy()
-
-
-def petsc2numpy_mat(petsc_mat):
-    """Allgather a PETSc.Mat."""
-    comm = petsc_mat.getComm()
-    local_mat = petsc_mat.getRedundantMatrix(
-        comm.size, PETSc.COMM_SELF)
-    return csr_matrix(
-        local_mat.getValuesCSR()[::-1],
-        shape=local_mat.getSize()
-    ).todense()
-
-
-@pytest.fixture
-def rng():
-    return RandomGenerator(PCG64(seed=13))
-
-
-@pytest.mark.skipcomplex
-@pytest.mark.parallel([1, 2])
-@pytest.mark.parametrize("degree", (1, 2), ids=["degree1", "degree2"])
-@pytest.mark.parametrize("dim", (0, 2, (2, 2)), ids=["scalar", "vec2", "tensor22"])
-@pytest.mark.parametrize("family", ("CG", "DG"))
-@pytest.mark.parametrize("mesh_type", ("interval", "square"))
-@pytest.mark.parametrize("backend_type", (PyOP2NoiseBackend, PetscNoiseBackend), ids=("pyop2", "petsc"))
-def test_white_noise(family, degree, mesh_type, dim, backend_type, rng):
-    """Test that white noise generator converges to a mass matrix covariance.
-    """
-
-    nx = 10
-    # Mesh dimension
-    if mesh_type == 'interval':
-        mesh = UnitIntervalMesh(nx)
-    elif mesh_type == 'square':
-        mesh = UnitSquareMesh(nx, nx)
-
-    # Variable rank
-    if not isinstance(dim, int):
-        V = TensorFunctionSpace(mesh, family, degree, shape=dim)
-    elif dim > 0:
-        V = VectorFunctionSpace(mesh, family, degree, dim=dim)
-    else:
-        V = FunctionSpace(mesh, family, degree)
-
-    # Finite element white noise has mass matrix covariance
-    M = inner(TrialFunction(V), TestFunction(V))*dx
-    covmat = petsc2numpy_mat(
-        assemble(M, mat_type='aij').petscmat)
-
-    generator = WhiteNoiseGenerator(
-        V, backend=backend_type(V, rng=rng))
-
-    # Test convergence as sample size increases
-    nsamples = [50, 100, 200, 400, 800]
-
-    samples = np.empty((V.dim(), nsamples[-1]))
-    for i in range(nsamples[-1]):
-        with generator.sample().dat.vec_ro as bv:
-            samples[:, i] = petsc2numpy_vec(bv)
-
-    covariances = [np.cov(samples[:, :ns]) for ns in nsamples]
-
-    # Measured covariance matrix should converge at a rate of sqrt(n).
-    # The number of samples is fairly small to keep the test cost
-    # lower so we only test the mean rate to a low tolerance.
-
-    errors = [np.linalg.norm(cov-covmat) for cov in covariances]
-    rate = -np.diff(np.log(errors))/np.diff(np.log(nsamples))
-
-    assert np.mean(rate) > 0.4
-
-
-@pytest.mark.skipcomplex
-@pytest.mark.parallel([1, 2])
-@pytest.mark.parametrize("dim", (0, 2, (2, 2)), ids=["scalar", "vec2", "tensor22"])
-@pytest.mark.parametrize("mesh_type", ("interval", "square"))
-def test_vom_white_noise(dim, mesh_type, rng):
-    """Test that white noise generator converges to a mass matrix covariance.
-    """
-
-    nx = 10
-    nv = 10
-    np.random.seed(13)
-    # Mesh dimension
-    if mesh_type == 'interval':
-        mesh = UnitIntervalMesh(nx)
-        points = np.random.random_sample((nv, 1))
-    elif mesh_type == 'square':
-        mesh = UnitSquareMesh(nx, nx)
-        points = np.random.random_sample((nv, 2))
-
-    vom = VertexOnlyMesh(mesh, points)
-
-    # Variable rank
-    if not isinstance(dim, int):
-        V = TensorFunctionSpace(vom, "DG", 0, shape=dim)
-    elif dim > 0:
-        V = VectorFunctionSpace(vom, "DG", 0, dim=dim)
-    else:
-        V = FunctionSpace(vom, "DG", 0)
-
-    # Finite element white noise has mass matrix covariance
-    M = inner(TrialFunction(V), TestFunction(V))*dx
-    covmat = petsc2numpy_mat(
-        assemble(M, mat_type='aij').petscmat)
-
-    backend = VOMNoiseBackend(V, rng)
-    generator = WhiteNoiseGenerator(V, backend=backend)
-
-    # Test convergence as sample size increases
-    nsamples = [50, 100, 200, 400, 800]
-
-    samples = np.empty((V.dim(), nsamples[-1]))
-    for i in range(nsamples[-1]):
-        with generator.sample().dat.vec_ro as bv:
-            samples[:, i] = petsc2numpy_vec(bv)
-
-    covariances = [np.cov(samples[:, :ns]) for ns in nsamples]
-
-    # Measured covariance matrix should converge at a rate of sqrt(n).
-    # The number of samples is fairly small to keep the test cost
-    # lower so we only test the mean rate to a low tolerance.
-
-    errors = [np.linalg.norm(cov-covmat) for cov in covariances]
-    rate = -np.diff(np.log(errors))/np.diff(np.log(nsamples))
-
-    assert np.mean(rate) > 0.4
-
-
-@pytest.mark.skipcomplex
-@pytest.mark.parallel([1, 2])
-@pytest.mark.parametrize("m", (0, 2, 4))
-@pytest.mark.parametrize("dim", (0, 2), ids=["scalar", "vector2"])
-@pytest.mark.parametrize("family", ("CG", "DG"))
-@pytest.mark.parametrize("mesh_type", ("interval", "square"))
-def test_covariance_inverse_action(m, family, mesh_type, dim):
-    """Test that covariance operator action and inverse are opposites.
-    """
-
-    nx = 20
-    if mesh_type == 'interval':
-        mesh = PeriodicUnitIntervalMesh(nx)
-        x, = SpatialCoordinate(mesh)
-        wexpr = cos(2*pi*x)
-    elif mesh_type == 'square':
-        mesh = PeriodicUnitSquareMesh(nx, nx)
-        x, y = SpatialCoordinate(mesh)
-        wexpr = cos(2*pi*x)*cos(4*pi*y)
-    elif mesh_type == 'cube':
-        mesh = PeriodicUnitCubeMesh(nx, nx, nx)
-        x, y, z = SpatialCoordinate(mesh)
-        wexpr = cos(2*pi*x)*cos(4*pi*y)*cos(pi*z)
-    if dim > 0:
-        V = VectorFunctionSpace(mesh, family, 1, dim=dim)
-        wexpr = as_vector([-1**(j+1)*wexpr for j in range(dim)])
-    else:
-        V = FunctionSpace(mesh, family, 1)
-
-    L = 0.1
-    sigma = 0.9
-
-    solver_parameters = {
-        'ksp_type': 'preonly',
-        'pc_type': 'lu',
-        'pc_factor_mat_solver_type': 'mumps'
-    }
-
-    form = 'IP' if family == 'DG' else 'CG'
-
-    B = AutoregressiveCovariance(
-        V, L, sigma, m, form=form,
-        solver_parameters=solver_parameters,
-        options_prefix="")
-
-    w = Function(V).project(wexpr)
-    wcheck = B.apply_action(B.apply_inverse(w))
-
-    tol = 1e-10
-
-    assert errornorm(w, wcheck) < tol
-
-
-@pytest.mark.skipcomplex
-@pytest.mark.parallel([1, 2])
-@pytest.mark.parametrize("m", (0, 2, 4))
-def test_covariance_inverse_action_hdiv(m):
-    """Test that covariance operator action and inverse are opposites
-    for hdiv spaces.
-    """
-
-    nx = 20
-    mesh = PeriodicUnitSquareMesh(nx, nx)
-    x, y = SpatialCoordinate(mesh)
-    wexpr = cos(2*pi*x)*cos(4*pi*x)
-
-    V = FunctionSpace(mesh, "BDM", 1)
-    wexpr = as_vector([-1**(j+1)*wexpr for j in range(2)])
-
-    L = 0.1
-    sigma = 0.9
-
-    solver_parameters = {
-        'ksp_type': 'preonly',
-        'pc_type': 'lu',
-        'pc_factor_mat_solver_type': 'mumps'
-    }
-
-    B = AutoregressiveCovariance(
-        V, L, sigma, m, form='IP',
-        solver_parameters=solver_parameters,
-        options_prefix="")
-
-    w = Function(V).project(wexpr)
-    wcheck = B.apply_action(B.apply_inverse(w))
-
-    tol = 1e-8
-
-    assert errornorm(w, wcheck) < tol
+@pytest.fixture(autouse=True)
+def autouse_set_test_tape(set_test_tape):
+    pass
 
 
 @pytest.mark.skipcomplex
@@ -269,133 +45,3 @@ def test_covariance_adjoint_norm(m, family):
     assert min(taylor['R0']['Rate']) > 0.95, taylor['R0']
     assert min(taylor['R1']['Rate']) > 1.95, taylor['R1']
     assert min(taylor['R2']['Rate']) > 2.95, taylor['R2']
-
-
-@pytest.mark.skipcomplex
-@pytest.mark.parallel([1, 2])
-@pytest.mark.parametrize("m", (0, 2, 4))
-@pytest.mark.parametrize("family", ("CG", "DG"))
-@pytest.mark.parametrize("operation", ("action", "inverse"))
-def test_covariance_mat(m, family, operation):
-    """Test that covariance mat and pc apply correct and opposite actions.
-    """
-    nx = 20
-    L = 0.2
-    sigma = 0.9
-
-    mesh = UnitIntervalMesh(nx)
-    coords, = SpatialCoordinate(mesh)
-
-    V = FunctionSpace(mesh, family, 1)
-
-    form = 'IP' if family == 'DG' else 'CG'
-
-    B = AutoregressiveCovariance(V, L, sigma, m, form=form)
-
-    mat = CovarianceMat(B, operation=operation)
-
-    expr = 2*pi*coords
-
-    if operation == 'action':
-        x = Function(V).project(expr).riesz_representation()
-        y = Function(V)
-        xcheck = x.copy(deepcopy=True)
-        ycheck = y.copy(deepcopy=True)
-
-        B.apply_action(xcheck, tensor=ycheck)
-
-    elif operation == 'inverse':
-        x = Function(V).project(expr)
-        y = Function(V.dual())
-        xcheck = x.copy(deepcopy=True)
-        ycheck = y.copy(deepcopy=True)
-
-        B.apply_inverse(xcheck, tensor=ycheck)
-
-    with x.dat.vec as xv, y.dat.vec as yv:
-        mat.mult(xv, yv)
-
-    # flip to primal space to calculate norms
-    if operation == 'inverse':
-        y = y.riesz_representation()
-        ycheck = ycheck.riesz_representation()
-
-    assert errornorm(ycheck, y)/norm(ycheck) < 1e-12
-
-    if operation == 'inverse':
-        y = y.riesz_representation()
-        ycheck = ycheck.riesz_representation()
-
-    ksp = PETSc.KSP().create()
-    ksp.setOperators(mat)
-
-    tol = 1e-8
-
-    petsctools.set_from_options(
-        ksp, options_prefix=str(operation),
-        parameters={
-            'ksp_monitor': None,
-            'ksp_type': 'richardson',
-            'ksp_max_it': 2,
-            'ksp_rtol': tol,
-            'pc_type': 'python',
-            'pc_python_type': 'firedrake.adjoint.CovariancePC',
-        }
-    )
-    x.zero()
-
-    with x.dat.vec as xv, y.dat.vec as yv:
-        with petsctools.inserted_options(ksp):
-            ksp.solve(yv, xv)
-
-    # CovarianceOperator operations should
-    # be exact inverses of each other.
-    assert ksp.its == 1
-
-    if operation == 'action':
-        x = x.riesz_representation()
-        xcheck = xcheck.riesz_representation()
-
-    assert errornorm(xcheck, x)/norm(xcheck) < 10*tol
-
-
-@pytest.mark.skipcomplex
-@pytest.mark.parametrize("family", ("CG", "DG"))
-def test_diffusion_form(family):
-    """Test that the provided diffusion forms converge to a known solution at the expected rate.
-    """
-    from firedrake.adjoint.covariance_operator import diffusion_form
-
-    def poisson_error(mesh, family):
-        V = FunctionSpace(mesh, family, 1)
-        x, y = SpatialCoordinate(mesh)
-
-        f = Function(V).interpolate((1+8*pi*pi)*cos(x*pi*2)*cos(y*pi*2))
-        uexact = cos(x*pi*2)*cos(y*pi*2)
-
-        nu = Constant(1)
-        u = TrialFunction(V)
-        v = TestFunction(V)
-
-        formulation = AutoregressiveCovariance.DiffusionForm(
-            "CG" if family == "CG" else "IP")
-
-        a = diffusion_form(u, v, nu, formulation)
-        L = inner(f, v)*dx
-
-        u = Function(V)
-        solve(a == L, u)
-
-        return errornorm(uexact, u)
-
-    base_nx = 16
-    nrefs = 3
-    base_mesh = UnitSquareMesh(base_nx, base_nx)
-    mh = MeshHierarchy(base_mesh, nrefs)
-
-    errors = [poisson_error(m, family) for m in mh]
-
-    # second order convergence
-    nxs = [base_nx*(2**n) for n in range(nrefs+1)]
-    rate = - np.diff(np.log(errors))/np.diff(np.log(nxs))
-    assert (rate > 1.9).all()