Integral conservation of adapted function?


133
views
0
4 months ago by
Hello,

it appears that adaptation of a function to a refined mesh does not preserve the integral of the function.

The following code produces two different results for me:

print(str(assemble(u0*dx)))

refined_mesh = refine(mesh, cell_markers)
adapt(u0,refined_mesh)
dx = Measure("dx",domain =refined_mesh)

print(str(assemble(u0*dx)))
​

Is this a typical behavior or am I doing something wrong?

Community: FEniCS Project
1
Not sure about the adapt method, but have you tried interpolating u0 onto the refined functionspace first? Also, providing a complete MWE would be appreciated.
written 4 months ago by pf4d  
Yes, I tried this, but it also failed.

I put together a 'MWE' based on the Cahn-Hilliard Demo.
The first few Iterations, the integral value stays the same and afterwards it changes, but only if I use Mesh adaptation.
import random
from dolfin import *

set_log_level(WARNING)

# Class representing the intial conditions
class InitialConditions(Expression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(mpi_comm_world()))
    def eval(self, values, x):
        values[0] = 0.63 + 0.02*(0.5 - random.random())
        values[1] = 0.0
    def value_shape(self):
        return (2,)

# Model parameters
lmbda  = 1.0e-02  # surface parameter
dt     = 5.0e-06  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson


# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

# Create mesh and build function space
mesh = UnitSquareMesh(96, 96)
mh = MeshHierarchy(mesh)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mh.finest(), P1*P1)

# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)
u.interpolate(u_init)
u0.interpolate(u_init)

# Compute the chemical potential df/dc
c = variable(c)
f    = 100*c**2*(1-c)**2
dfdc = diff(f, c)

# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu

# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Create nonlinear problem and Newton solver
problem = NonlinearVariationalProblem(L, u, J=a)
solver =  NonlinearVariationalSolver(problem)

# Output file
file = File("output.pvd", "compressed")

# Step in time
t = 0.0
T = 50*dt
while (t < T):
    t += dt
    u0.vector()[:] = u.vector()

    mh = MeshHierarchy(mesh)
    #to access values
    c0, mu0 = u0.split()

    for i in range(2):
        c1 = project(c0,FunctionSpace(mh.finest(),P1))
        cell_markers_refine = CellFunction("bool", mh.finest())
        cell_markers_refine.set_all(False)
        for cell in  cells(mh.finest()):
            point = cell.midpoint()
            tmp = abs(c1([point.x(), point.y()]))
            if (cell.h() > .01) and (tmp > 0.1 and tmp < 0.9):
                cell_markers_refine[cell] = True

        mh = mh.refine(cell_markers_refine)

        #adapt to refined mesh
        dx = Measure("dx",domain = mh.finest())
        ds = Measure("ds",domain = mh.finest())

        adapt(u0,mh.finest())
        adapt(u,mh.finest())
        adapt(problem, mh.finest())
        u0        = u0.child()
        u         = u.child()
        problem   = problem.child()
        solver    = NonlinearVariationalSolver(problem.leaf_node())


    solver.solve()

    c0, mu0 = u0.split()
    print(str(assemble(c0*dx)))
    # file << (u.split()[0], t)
​
written 3 months ago by Lukas O.  

1 Answer


0
5 weeks ago by
To answer my own question: Since my plan was to use "adapt" really to adapt the mesh, i.e. not only refine it, but also achieve some kind of coarsening,  I started in every timestep from the coarsest mesh and refined it. Afterwards I adapted the solution to the refined mesh. However, that means that it can occur that the solution needs to be interpolated from finer to coarser parts of the grid. Then there is an interpolation error leading to the loss of integral conservation.

As a follow up question: Is it possible to prevent this, meaning using the integral conservation as constraint for adapting the function?
Please login to add an answer/comment or follow this question.

Similar posts:
Search »