Galerkin method no longer resolving the sharp gradient layer

General
1

Hello dr @dokken I am having some issues with dolfinx. I had previously done some work on one of my problems and the galerkin method was completely resolving the layer (images attached) for the same mesh size and diffusion coefficient epsilon. it is no longer resolving the fine layer of the solution. I am attaching my previous images I had and new ones as well. I am supposed to do some error analysis now but the I am not getting what I was expecting. I am not sure whether it is a bug in new 0.8 version or something else.
Code:

import numpy as np

from dolfinx.fem import (Constant, Function, functionspace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx import fem, io, mesh
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, dx, ds, grad, inner
from petsc4py.PETSc import ScalarType


def solve_adv_diff(N=100, eps= ScalarType(1), degree=1):

    # define mesh and function space
    msh = create_unit_square(MPI.COMM_WORLD,N,N)
    V = functionspace(msh, ("Lagrange", degree))
    u = TrialFunction(V)
    v = TestFunction(V)
    x = SpatialCoordinate(msh)

    # define boundary conditions
    # bottom left_facet
    facet_1 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.logical_and(np.isclose(x[1], 0),x[0]<=0.5))
    dof_1 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_1)
    bc_1 = dirichletbc(ScalarType(1), dof_1,V)

    # left facet
    facet_2 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.isclose(x[0], 0))
    dof_2 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_2)
    bc_2 = dirichletbc(ScalarType(1), dof_2,V)

    # bottom right_facet
    facet_3 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.logical_and(np.isclose(x[1], 0),x[0]>0.5))
    dof_3 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_3)
    bc_3 = dirichletbc(ScalarType(0), dof_3,V)
    bcs = [bc_1,bc_2, bc_3]

    # define variational form
    vec_f = Constant(msh,ScalarType((np.cos(np.pi/3),np.sin(np.pi/3))))
    a = eps*dot(grad(u), grad(v))*dx + v*dot(vec_f,grad(u))*dx
    f = Constant(msh, ScalarType(0))
    L = f * v * dx

    # solving problem
    problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
    return problem.solve() 

u_galerkin = solve_adv_diff(100,1)
# writing mesh and solution for visualization in paraview
from dolfinx import io
with io.XDMFFile(u_galerkin.function_space.mesh.comm, "u_h100_eps1.xdmf", "w") as file:
    file.write_mesh(u_galerkin.function_space.mesh)
    file.write_function(u_galerkin)

h_100eps1
h_100eps11304Ă—640 34.5 KB
(New image)
100
1001304Ă—640 36.7 KB
(old image)

2

@dokken is a wizard :mage:, but if you don’t provide him the full code it’s unlikely that he (or anyone else) will be able to help :stuck_out_tongue:

Looking at your pictures, I would bet that either you are using a much smaller value for N, or you are using a much larger value of eps.

3

@francesco-ballarin I have provided the full code for getting the solution. Old pic was generated with the same mesh size and even smaller eps it still resolved the fine layer. eps is 1 here. I even used N = 1000 still didn’t resolve that’s why I am asking

4

That’s not the full code. Where are the imports? How are we supposed to know if you are calling solve_adv_diff with the default arguments or with other arguments?

6 
def error_L2(uh, u_ex, degree_raise = 3):
    # Create higher order function space
    tol = 1e-6
    degree = u_ex.function_space.ufl_element().degree
    family = "CG"
    mesh = u_ex.function_space.mesh
    W = functionspace(mesh, (family, degree + degree_raise))
    u_fine = Function(W)
    u_fine.interpolate(u_ex) 
    u_coarse = Function(W)
    # Interpolate approx. solution in fine solution function space

    if uh.function_space.mesh != mesh:
        from dolfinx.fem import create_nonmatching_meshes_interpolation_data
        interpolation_data = create_nonmatching_meshes_interpolation_data(
            u_coarse.function_space.mesh._cpp_object,
            u_coarse.function_space.element,
            uh.function_space.mesh._cpp_object,tol)
        u_coarse.interpolate(uh, nmm_interpolation_data=interpolation_data)
    else:
        u_coarse.interpolate(uh)

    # Compute the error in fine solution function space
    e_W = Function(W)
    e_W.x.array[:] = u_fine.x.array - u_coarse.x.array

    # Integrate the error
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)


num_elements_list = [25, 50,100]
eps_list = [0.001,0.01,0.1,1]

print("Compare to finely resolved approximation")
for num_elements in num_elements_list:
  for eps in eps_list:
    u_approx = solve_adv_diff(num_elements,eps)
    u_fine = solve_adv_diff(200,eps)
    print(f"{num_elements}: {eps}: {error_L2(u_approx, u_fine)}")

error output goes like this:

Compare to finely resolved approximation
25: 0.001: 0.013683519989143383
25: 0.01: 0.005832569416414266
25: 0.1: 0.004613618901075035
25: 1: 0.0044905001329669635
50: 0.001: 0.02203068334229291
50: 0.01: 0.013594735336098614
50: 0.1: 0.010021462565642979
50: 1: 0.008782639246786884
100: 0.001: 0.007560060292445784
100: 0.01: 0.004730443243876182
100: 0.1: 0.003561693981855338
100: 1: 0.0031668502681140758

example: we can see error increases when we increase N from 25 to 50 for eps = 1

7

You are really making it hard for people to help you. You should have posted from the very beginning (or, when asked for the first time) the full code with imports, and not spread bits and pieces in three different posts. The error computation is unnecessary.

Regardless, you have

u_galerkin = solve_adv_diff(100,1)

i.e. you have a problem with diffusivity equal to 1 and a transport field of magnitude 1. Given those values, the first picture is what I would expect from a physical point of view, so there is nothing wrong with the code. Go back to your old code and check the value of diffusivity that you used there, because I am confident that it wasn’t one (or the transport field had a different magnitude).

8

@francesco-ballarin I am really sorry for that. I have edited it again and posted the error computation we can see error increasing when we increase N from 25 to 50 for eps = 1. When I computed previously I used eps around 0.01 with N = 100 and it was still resolving with some little to no oscillations. I don’t have the previous code now. I have updated it to what I posted. Can you please comment on why error is increasing?? Am I doing anything wrong with error computation.

9

If I were you, I would avoid comparing the error obtained with an odd value of N to the error with an even value. The case of odd N is “dangerous” because the point (0.5, 0) is in the middle of a facet, and so are you are leaving to an internal implementation detail in dolfinx whether that facet will be marked as associated to BC equal to zero or BC equal to 1. If N is even, the point (0.5, 0) will surely be a vertex of a boundary facet.

In the code below I only compare even cases. The errors do decrease with refinement, and I’ll leave it to you to verify if they do so with the expected convergence rate (and, if the don’t, why).

Note that I switched the loops in the error computation, so that for a given value of eps you only have to compute the fine solution once.

import numpy as np

from dolfinx.fem import (Constant, Function, functionspace,
                         assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx import fem, io, mesh
from dolfinx.mesh import create_unit_square
from mpi4py import MPI
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, dx, ds, grad, inner
from petsc4py.PETSc import ScalarType


def solve_adv_diff(N=100, eps= ScalarType(1), degree=1):

    # define mesh and function space
    msh = create_unit_square(MPI.COMM_WORLD,N,N)
    V = functionspace(msh, ("Lagrange", degree))
    u = TrialFunction(V)
    v = TestFunction(V)
    x = SpatialCoordinate(msh)

    # define boundary conditions
    # bottom left_facet
    facet_1 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.logical_and(np.isclose(x[1], 0),x[0]<=0.5))
    dof_1 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_1)
    bc_1 = dirichletbc(ScalarType(1), dof_1,V)

    # left facet
    facet_2 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.isclose(x[0], 0))
    dof_2 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_2)
    bc_2 = dirichletbc(ScalarType(1), dof_2,V)

    # bottom right_facet
    facet_3 = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                        marker=lambda x: np.logical_and(np.isclose(x[1], 0),x[0]>0.5))
    dof_3 = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facet_3)
    bc_3 = dirichletbc(ScalarType(0), dof_3,V)
    bcs = [bc_1,bc_2, bc_3]

    # define variational form
    vec_f = Constant(msh,ScalarType((np.cos(np.pi/3),np.sin(np.pi/3))))
    a = eps*dot(grad(u), grad(v))*dx + v*dot(vec_f,grad(u))*dx
    f = Constant(msh, ScalarType(0))
    L = f * v * dx

    # solving problem
    problem = LinearProblem(a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
    return problem.solve()

def error_L2(uh, u_ex, degree_raise = 3):
    # Create higher order function space
    tol = 1e-16
    degree = u_ex.function_space.ufl_element().degree
    family = "CG"
    mesh = u_ex.function_space.mesh
    W = functionspace(mesh, (family, degree + degree_raise))
    u_fine = Function(W)
    u_fine.interpolate(u_ex)
    u_coarse = Function(W)
    # Interpolate approx. solution in fine solution function space

    if uh.function_space.mesh != mesh:
        from dolfinx.fem import create_nonmatching_meshes_interpolation_data
        interpolation_data = create_nonmatching_meshes_interpolation_data(
            u_coarse.function_space.mesh._cpp_object,
            u_coarse.function_space.element,
            uh.function_space.mesh._cpp_object,
            tol)
        u_coarse.interpolate(uh, nmm_interpolation_data=interpolation_data)
    else:
        u_coarse.interpolate(uh)

    # Compute the error in fine solution function space
    e_W = Function(W)
    e_W.x.array[:] = u_fine.x.array - u_coarse.x.array

    # Integrate the error
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)


num_elements_list = [24, 50, 100]
# num_elements_list = [l + 1 for l in num_elements_list]  # DANGEROUS!!
eps_list = [0.001,0.01,0.1,1]

print("Compare to finely resolved approximation")
for eps in eps_list:
    u_fine = solve_adv_diff(200,eps)
    for num_elements in num_elements_list:
        u_approx = solve_adv_diff(num_elements,eps)
        print(f"{num_elements}: {eps}: {error_L2(u_approx, u_fine)}")
1 Like
10

I appreciate your patience with me. I will surely check the convergence rate as well. Thanks for pointing out this odd and even issue. I am truly grateful!