Different results on both sides of a variational problem

variational formulation
1

Hi everyone,
I am trying to solve a non-linear temporal problem with two unknowns (eta and c). The domain is one-dimensional and boundary conditions and a given initial condition apply. The strong formulation of the problem is the following:

eq1
eq1788×128 8.15 KB

and the weak formulation is the following solving with an implicit temporal scheme.

eq2
eq2709×178 8.7 KB

When I introduce the values of the unknowns in the strong formulation, shouldn’t the equation be satisfied? However the results of the two parts of the equation are different.

im1
im1504×744 56.9 KB

This is tehe code:

import dolfinx.generation
import dolfinx.fem
import dolfinx.mesh
import dolfinx.nls
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from matplotlib import pyplot as plt

#Mesh creation
nx = 1000     #number of nodes
mesh = dolfinx.generation.IntervalMesh(MPI.COMM_WORLD, nx-1 ,[0,150e-6])

#Mixed function space creation
P1 = ufl.FiniteElement("CG", mesh.ufl_cell(), 2)
ME = dolfinx.fem.FunctionSpace(mesh,ufl.MixedElement([P1, P1]))

#Constant definition
D1 = dolfinx.fem.Constant(mesh, 8.5e-10)  
A = dolfinx.fem.Constant(mesh, 5.35e7)    
L = dolfinx.fem.Constant(mesh, 2)         
l = dolfinx.fem.Constant(mesh, 5e-6)      
sigma = dolfinx.fem.Constant(mesh, 10)    
alpha = dolfinx.fem.Constant(mesh, 2.94)  
c_le = dolfinx.fem.Constant(mesh, 0.03566433566433566)
c_se = dolfinx.fem.Constant(mesh, 1)
alpha_eta = dolfinx.fem.Constant(mesh, 3.006406382595866e-06)
w_ = dolfinx.fem.Constant(mesh,2078893.9366884497)

#Function definition
def H(eta):
    return -2.*eta**3 + 3.*eta**2
def derH(eta):
    return 6.*(eta-eta**2)
def G(eta):
    return (1.-eta)**2*eta**2
def derG(eta):
    return 2*eta-6*eta**2+4*eta**3

#Boundary conditions
uD = dolfinx.fem.Function(ME)
uD.sub(0).interpolate(lambda x: 1 - x[0]/150e-6)
uD.sub(1).interpolate(lambda x: 1 - x[0]/150e-6)
uD.x.scatter_forward()
boundary_facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim - 1, lambda x: np.full(x.shape[1], True, dtype=bool))
boundary_dofs = dolfinx.fem.locate_dofs_topological(ME, mesh.topology.dim - 1, boundary_facets)
bc = dolfinx.fem.DirichletBC(uD, boundary_dofs)

# Define test functions
q, v = ufl.TestFunctions(ME)
# Define functions
u = dolfinx.fem.Function(ME)
u0 = dolfinx.fem.Function(ME)
# Split mixed functions
c, eta = ufl.split(u)
c0, eta0 = ufl.split(u0)

# Define initial condition
u.sub(0).interpolate(lambda x: 1/(1+ufl.e**(2e6*(x[0]-0.90*150e-6))))
u.sub(1).interpolate(lambda x: 1/(1+ufl.e**(1.6e6*(x[0]-0.90*150e-6))))
u.x.scatter_forward()
u0.sub(0).interpolate(lambda x: 1/(1+ufl.e**(2e6*(x[0]-0.90*150e-6))))
u0.sub(1).interpolate(lambda x: 1/(1+ufl.e**(1.6e6*(x[0]-0.90*150e-6))))
u0.x.scatter_forward()

# Define temporal parameters
t = 0 # Start time
T = 10 # Final time
dt = 0.5 # Interval time
num_steps = int(T/dt)
k = dolfinx.fem.Constant(mesh, dt)

# Variational problem
F0 = ufl.inner(eta, v)*ufl.dx - ufl.inner(eta0, v)*ufl.dx + k*2*L*A*(c_le-1)*ufl.inner((c+c_le*(H(eta)-1.)-H(eta))*derH(eta), v)*ufl.dx + k*L*alpha_eta*ufl.inner(ufl.grad(eta), ufl.grad(v))*ufl.dx  + k*L*w_*ufl.inner(derG(eta),v)*ufl.dx
F1 = ufl.inner(c,q)*ufl.dx - ufl.inner(c0,q)*ufl.dx + k*ufl.inner(D1*ufl.grad(c), ufl.grad(q))*ufl.dx - k*(1-c_le)*derH(eta)*ufl.inner(D1*ufl.grad(eta), ufl.grad(q))*ufl.dx
F = F0 + F1
problem = dolfinx.fem.NonlinearProblem(F, u, bcs=[bc])
solver = dolfinx.nls.NewtonSolver(MPI.COMM_WORLD, problem)
solver.atol = 1e-12
solver.rtol = 1e-12
solver.convergence_criterion = "incremental"
solver.report = True
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

tol = 1e-8
while (t < (T-tol)):
    t += dt
    r = solver.solve(u)
    u.x.scatter_forward()
    if t < (T-tol): #To keep u0 in last step
        u.vector.copy(result=u0.vector)

#Project both sides of the equation to compare       
V = dolfinx.fem.FunctionSpace(mesh, ("CG", 2))

def project(expression, V):
    u_, v_ = ufl.TrialFunction(V), ufl.TestFunction(V)
    a_p = ufl.inner(u_, v_) * ufl.dx
    L_p = ufl.inner(expression, v_) * ufl.dx
    projection = dolfinx.fem.LinearProblem(a_p, L_p) #dolfinx.fem.LinearProblem(a, L, [bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"}) #dolfinx.fem.LinearProblem(a, L, bcs=[], u=None, petsc_options={},form_compiler_parameters={"ksp_type": "preonly", "pc_type": "lu"},jit_parameters={})
    return projection.solve()

fb1 = - k*derH(eta)*(c-H(eta)*(c_se-c_le)-c_le)*2*A*L*(c_le-c_se) - k*L*w_*derG(eta) + k*L*alpha_eta*ufl.div(ufl.grad(eta))

print('First term integral:',dolfinx.fem.assemble_scalar((eta-eta0)* ufl.dx))
print('Second term integral:',dolfinx.fem.assemble_scalar(fb1* ufl.dx))

projection1 = project(eta-eta0,V)#projection of first term
projection2 = project(fb1,V)#projection of second term

plt.figure()
plt.plot(mesh.geometry.x[:,0], projection1.compute_point_values().real[:,0])
plt.title('First term')
plt.xlabel('x[m]')
plt.ylabel('value[-]')
plt.grid()

plt.figure()
plt.plot(mesh.geometry.x[:,0], projection2.compute_point_values().real[:,0])
plt.title('Second term')
plt.xlabel('x[m]')
plt.ylabel('value[-]')
plt.grid()

I would be very grateful if anyone knows why this is happening. :grin:

2

When I introduce the values of the unknowns in the strong formulation, shouldn’t the equation be satisfied?

You don’t expect the strong form of the PDE to be satisfied exactly by the finite element solution. Considering the weighted-residual derivation, you might initially expect that the strong problem would be satisfied in an average sense over the domain, which is what you’re testing in your MWE. However, global constants are not part of the test function space, due to the DirichletBC, and the form fb1 is missing distributional contributions to the finite element solution’s second derivatives (as I discussed in more detail in my answer here).

1 Like
3

First of all thank you very much for the great answer, it has helped me to clarify some concepts :grinning:

I would like to ask you a couple of questions:

  1. On the one hand I thought that FEniCS/FEniCSX used the Galerking method to obtain the algebraic equations to be solved. Is the Galerking method an example of the method of weighted residues? What does FEniCS use to obtain the algebraic equations?

  2. If I have understood correctly, one of the problems is that I am trying to project the Laplacian of one of the unknowns and when I have to evaluate it in the interior of each element this value is zero if the elements that are piecewise linear and in other elements the Dirac delta would be missing. I have tried to increase the degree of the elements and each time the results are more similar, obtaining the following result with elements of degree 4:

Imagen13
Imagen13424×621 60 KB

Would there be any other way to correctly evaluate the values of the Laplacian term that the solver has taken into account when giving a time step?

4

On the one hand I thought that FEniCS/FEniCSX used the Galerking method to obtain the algebraic equations to be solved. Is the Galerking method an example of the method of weighted residues? What does FEniCS use to obtain the algebraic equations?

FEniCS uses whatever variational method you choose to implement! It looks like your MWE is using the standard Galerkin method (sometimes referred to as “Bubnov–Galerkin” in settings where “Galerkin” refers to variational methods generally), but you could also implement various Petrov–Galerkin, discontinuous Galerkin, weakly-consistent, or other methods.

Would there be any other way to correctly evaluate the values of the Laplacian term that the solver has taken into account when giving a time step?

Extracting higher derivatives from low-order finite element solutions can be tricky. You might take a look at the discussion here for some possibilities.

1 Like