adaptive_example_w_mixedelement.py

import dolfin as dolf
import numpy as np

LOG_LEVEL = 30
dolf.set_log_level(LOG_LEVEL)

mesh = dolf.UnitIntervalMesh(100)
n = dolf.FacetNormal(mesh)
wall_left = "near(x[0], 0)"
wall_right = "near(x[0], 1.)"
all_bound = "on_boundary"

P2 = dolf.FiniteElement("Lagrange", mesh.ufl_cell(), 2)
P1 = dolf.FiniteElement("Lagrange", mesh.ufl_cell(), 2)
P_mixed = dolf.MixedElement([P2, P1])

W = dolf.FunctionSpace(mesh, P_mixed)
w_trial = dolf.TrialFunction(W)
w_test = dolf.TestFunction(W)

u, p = dolf.split(w_trial)
v, q = dolf.split(w_test)

dt = 1.0e-3
NSteps = int(31.0e2)


def B_Form(trial_func, test_func):
    u, p = trial_func
    v, q = test_func
    return (p * q + u * v) * dolf.dx


def A_Form_viscous(trial_func, test_func):
    u, p = trial_func
    v, q = test_func
    mu = dolf.Constant(1.0e-2)
    return (-u.dx(0) * q - (-p + mu * u.dx(0)) * v.dx(0)) * dolf.dx


def unsteady_solve(control_velocity=None, plotting=False):

    state_0 = dolf.Function(W)
    (u0, p0) = state_0.split(True)

    w = dolf.Function(W)

    for i in range(NSteps):

        a = B_Form((u, p), (v, q)) - dolf.Constant(dt * 0.5) * A_Form_viscous(
            (u, p), (v, q)
        )

        L = B_Form((u0, p0), (v, q)) + dolf.Constant(dt * 0.5) * A_Form_viscous(
            (u0, p0), (v, q)
        )
        L -= dolf.Constant(dt * 1.0) * n[0] * v * dolf.ds

        control_right = dolf.DirichletBC(W.sub(0), dolf.Constant(0.0), wall_right)

        bcs = [control_right]

        tol = 1.0e-5
        # How should I specify a goal here?
        # M = w * dolf.dx() ?
        # My guess was:
        (_u, _p) = w.split()
        M = (_u**2. + _p**2)*dolf.dx()

        problem = dolf.LinearVariationalProblem(a, L, w, bcs)
        solver = dolf.AdaptiveLinearVariationalSolver(problem, M)
        solver.parameters["error_control"]["dual_variational_solver"][
            "linear_solver"
        ] = "cg"
        solver.parameters["error_control"]["dual_variational_solver"][
            "symmetric"
        ] = True
        solver.solve(tol)
        solver.summary()

        (u0, p0) = w.split()


unsteady_solve()