Refined Error

[zhyxue0]

Hello, everyone.
I want to use the smallest refined solution as a Reference solution to calculate the error. But the program reported an error（Whether running in parallel or on a single core）.

Here is my code:

from firedrake import *
import math 
Test_Type = "Temporal_Convergence"


FinalTime = 0.5    
num_refs = 4
u_solution = []
for idx in range(num_refs):
    Lx = 2*math.pi; 
    Ly = 2*math.pi
    if Test_Type == "Spatial_Convergence":
        Nx = 4*2**(idx)
        N_step = 1024
    elif Test_Type == "Temporal_Convergence":
        Nx = 64
        N_step = 2**(idx)
        
    Ny = Nx
    Hx = Lx/Nx
    
    
    dt = FinalTime / N_step
    dtc = Constant(dt)
    
    t = 0
    tc = Constant(t)
    
    mesh = PeriodicRectangleMesh(Nx, Ny, Lx, Ly)
    x, y  = SpatialCoordinate(mesh)

    V = FunctionSpace(mesh, "Lagrange", 1)
    u_trial = TrialFunction(V)
    v_test = TestFunction(V)
    
    u_exact =  sin(x)*sin(y)*exp(-2*tc)
    u = Function(V).interpolate(u_exact)
    
    LHS = (
           inner(u_trial, v_test)*dx
          -0.5*dtc* inner(grad(u_trial), grad(v_test)) * dx 
           )
    RHS = (
           inner(u, v_test)*dx
          +0.5*dtc* inner(grad(u), grad(v_test)) * dx 
          )
    
    
    Solution_space = Function(V)
    
    solver_parameters={
                    "ksp_type": "preonly",
                    "pc_type": "lu",
                    "pc_factor_mat_solver_type": "mumps",
                    }    
    problem = LinearVariationalProblem(LHS, RHS, Solution_space)
    solver = LinearVariationalSolver(problem, solver_parameters= solver_parameters)

    Step = 0 
    for Step in range(N_step):
        t += dt
        tc.assign(t)
        solver.solve()       
        u.assign(Solution_space)

    u_solution.append(u)
   
# Use the smallest refined solution as a Reference solution to calculate the error
u_e_ref_L2 = []    
for i, u in enumerate(u_solution[:-1]):
    u_e_ref_L2.append(norm(u - u_solution[-1]), "L2")
    PETSc.Sys.Print("u_e_ref_L2={:}",u_e_ref_L2)

Here is the report:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[21], line 72
     67 u_e_ref_L2 = []    
     68 for i, u in enumerate(u_solution[:-1]):
     69     # u_ref = u_solution[i+1]
     70     # u_inter = project(u, u_ref.function_space())
     71     # error = errornorm(u_ref, u_inter)
---> 72     u_e_ref_L2.append(norm(u - u_solution[-1]), "L2")
     73     PETSc.Sys.Print("u_e_ref_L2={:}",u_e_ref_L2)

File petsc4py/PETSc/Log.pyx:115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py/PETSc/Log.pyx:116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File /home/firedrake/firedrake/src/firedrake/firedrake/norms.py:80, in norm(v, norm_type, mesh)
     77 else:
     78     raise RuntimeError("Unknown norm type '%s'" % norm_type)
---> 80 return assemble((expr**(p/2))*dx)**(1/p)

File /home/firedrake/firedrake/src/ufl/ufl/measure.py:408, in Measure.__rmul__(self, integrand)
    406 domain = self.ufl_domain()
    407 if domain is None:
--> 408     domains = extract_domains(integrand)
    409     if len(domains) == 1:
    410         domain, = domains

File /home/firedrake/firedrake/src/ufl/ufl/domain.py:340, in extract_domains(expr)
    338 for t in traverse_unique_terminals(expr):
    339     domainlist.extend(t.ufl_domains())
--> 340 return sorted(join_domains(domainlist))

TypeError: '<' not supported between instances of 'MeshGeometry' and 'MeshGeometry'

[ksagiyam]

In this context Firedrake only works with a single domain, so you will first need to project the solution on the coarse mesh onto the finest mesh, and then compare that projected solution with u_solution[-1].

[zhyxue0]

In this context Firedrake only works with a single domain, so you will first need to project the solution on the coarse mesh onto the finest mesh, and then compare that projected solution with u_solution[-1].

This is the heat equation. My purpose is to test the accuracy of time, so the mesh is not refined. What is refined is the time step (as you can see, I fixed the space).

if Test_Type == "Spatial_Convergence":
        Nx = 4*2**(idx)
        N_step = 1024
elif Test_Type == "Temporal_Convergence":
        Nx = 64
        N_step = 2**(idx)

[connorjward]

I think the issue may be because each of these functions is created on a different mesh. Could you create the mesh outside the loop and reuse it?

[zhyxue0]

I think the issue may be because each of these functions is created on a different mesh. Could you create the mesh outside the loop and reuse it?

Thanks. It works. But if I want to test spatial accuracy, can I only project the solution onto the finest mesh? But it will introduce projection errors. How to solve it ?

[colinjcotter]

Thanks. It works. But if I want to test spatial accuracy, can I only project the solution onto the finest mesh? But it will introduce projection errors. How to solve it ? For most finite element spaces on nested meshes, the coarse mesh space is contained within the fine space. Then projection results in the same function but now expanded on the basis of the fine space mesh. So there is no error.

[zhyxue0]

Thanks. It works. But if I want to test spatial accuracy, can I only project the solution onto the finest mesh? But it will introduce projection errors. How to solve it ? For most finite element spaces on nested meshes, the coarse mesh space is contained within the fine space. Then projection results in the same function but now expanded on the basis of the fine space mesh. So there is no error.

Thanks. @colinjcotter. But the BDM element seems unable to do the projection operations?
Here is my code:

from firedrake import *
import math 
from firedrake.petsc import PETSc  # for PETSc.Sys.Print
import numpy as np
# Test_Type = "Temporal_Convergence"
Test_Type = 'Spatial_Convergence'

FinalTime = 0.5    
num_refs = 2

u_solution = []
mesh_ref = []
u_e_exact = []


Lx = 2*math.pi; 
Ly = 2*math.pi;

for idx in range(num_refs):
    # Lx = 2*math.pi; 
    # Ly = 2*math.pi
    if Test_Type == "Spatial_Convergence":
        Nx = 4*2**(idx)
        N_step = 1024
    elif Test_Type == "Temporal_Convergence":
        # Nx = 64
        N_step = 2**(idx)
        
    Ny = Nx
    Hx = Lx/Nx
    
    
    mesh = PeriodicRectangleMesh(Nx, Ny, Lx, Ly)
    x, y  = SpatialCoordinate(mesh)
    
    dt = FinalTime / N_step
    dtc = Constant(dt)
    
    t = 0
    tc = Constant(t)
    


    V = FunctionSpace(mesh, "BDM", 1)
    u_trial = TrialFunction(V)
    v_test = TestFunction(V)
    
    u_exact =  as_vector(
                        [cos(2*x)*sin(2*y)*exp(-8*tc), 
                         cos(2*x)*sin(2*y)*exp(-8*tc)]
                        )
    
    u = Function(V).interpolate(u_exact)
    
    LHS = (
           inner(u_trial, v_test)*dx
           + 0.5*dtc* inner(grad(u_trial), grad(v_test)) * dx 
           )
    
    RHS = (
           inner(u, v_test)*dx
          -0.5*dtc* inner(grad(u), grad(v_test)) * dx 
          )
    
    
    Solution = Function(V)
    
    solver_parameters={
                    "ksp_type": "preonly",
                    "pc_type": "lu",
                    "pc_factor_mat_solver_type": "mumps",
                    }    
    problem = LinearVariationalProblem(LHS, RHS, Solution)
    solver = LinearVariationalSolver(problem, solver_parameters= solver_parameters)

    Step = 0 
    for Step in range(N_step):
        t += dt
        solver.solve()       
        u.assign(Solution)
        # print(Step, norm(u-u_exact, "L2"))
    tc.assign(FinalTime)
    u_e_exact.append(norm(u-u_exact, "L2"))
    
    u_solution.append(u)
    
    
# Use the smallest refined solution as a Reference solution to calculate the error
u_e_ref_L2 = []    
u_ref = u_solution[-1]
for i, u in enumerate(u_solution[:-1]):
    u_inter = project(u, u_ref.function_space())
    u_e_ref_L2.append(norm(u_inter - u_solution[-1], "L2"))
    # PETSc.Sys.Print("u_e_ref_L2={:}",u_e_ref_L2)
    # PETSc.Sys.Print("u_e_ref_L2={:}",u_e_ref_L2)

Error = np.array(u_e_ref_L2)
Order_ref = np.log(Error[:-1]/Error[1:]) / np.log(2)
u_e_exact_np = np.array(u_e_exact)
Order_exact = np.log(u_e_exact_np[:-1]/u_e_exact_np[1:]) / np.log(2)

Here is the report:

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
Cell In[10], line 99
     97 u_ref = u_solution[-1]
     98 for i, u in enumerate(u_solution[:-1]):
---> 99     u_inter = project(u, u_ref.function_space())
    100     u_e_ref_L2.append(norm(u_inter - u_solution[-1], "L2"))
    101     # PETSc.Sys.Print("u_e_ref_L2={:}",u_e_ref_L2)
    102     # PETSc.Sys.Print("u_e_ref_L2={:}",u_e_ref_L2)

File petsc4py/PETSc/Log.pyx:115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py/PETSc/Log.pyx:116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File /home/firedrake/firedrake/src/firedrake/firedrake/adjoint/projection.py:34, in annotate_project..wrapper(*args, **kwargs)
     31             block = ProjectBlock(args[0], V, output, bcs, ad_block_tag=ad_block_tag, **sb_kwargs)
     33 with stop_annotating():
---> 34     output = project(*args, **kwargs)
     36 if annotate:
     37     tape = get_working_tape()

File /home/firedrake/firedrake/src/firedrake/firedrake/projection.py:78, in project(v, V, bcs, solver_parameters, form_compiler_parameters, use_slate_for_inverse, name, ad_block_tag)
     48 @PETSc.Log.EventDecorator()
     49 @annotate_project
     50 def project(v, V, bcs=None,
   (...)
     54             name=None,
     55             ad_block_tag=None):
     56     """Project a UFL expression into a :class:`.FunctionSpace`
     57     It is possible to project onto the trace space 'DGT', but not onto
     58     other trace spaces e.g. into the restriction of CG onto the facets.
   (...)
     73     then ``v`` is projected into a new :class:`.Function` and that
     74     :class:`.Function` is returned."""
     76     val = Projector(v, V, bcs=bcs, solver_parameters=solver_parameters,
     77                     form_compiler_parameters=form_compiler_parameters,
---> 78                     use_slate_for_inverse=use_slate_for_inverse).project()
     79     val.rename(name)
     80     return val

File /home/firedrake/firedrake/src/firedrake/firedrake/projection.py:181, in ProjectorBase.project(self)
    180 def project(self):
--> 181     self.apply_massinv(self.target, self.rhs)
    182     return self.target

File /home/firedrake/firedrake/src/firedrake/firedrake/projection.py:240, in SupermeshProjector.rhs(self)
    237 @property
    238 def rhs(self):
    239     with self.source.dat.vec_ro as u, self.residual.dat.vec_wo as v:
--> 240         self.mixed_mass.mult(u, v)
    241     return self.residual

File /home/firedrake/firedrake/src/PyOP2/pyop2/utils.py:62, in cached_property.__get__(self, obj, cls)
     60 if obj is None:
     61     return self
---> 62 obj.__dict__[self.__name__] = result = self.fget(obj)
     63 return result

File /home/firedrake/firedrake/src/firedrake/firedrake/projection.py:234, in SupermeshProjector.mixed_mass(self)
    231 @cached_property
    232 def mixed_mass(self):
    233     from firedrake.supermeshing import assemble_mixed_mass_matrix
--> 234     return assemble_mixed_mass_matrix(self.source.function_space(),
    235                                       self.target.function_space())

File petsc4py/PETSc/Log.pyx:115, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File petsc4py/PETSc/Log.pyx:116, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func()

File /home/firedrake/firedrake/src/firedrake/firedrake/supermeshing.py:76, in assemble_mixed_mass_matrix(V_A, V_B)
     70     if V_A.ufl_element().mapping() != "identity" or V_B.ufl_element().mapping() != "identity":
     71         msg = """
     72 Sorry, only implemented for affine maps for now. To do non-affine, we'd need to
     73 import much more of the assembly engine of UFL/TSFC/etc to do the assembly on
     74 each supermesh cell.
     75 """
---> 76         raise NotImplementedError(msg)
     78     mesh_A = V_A.mesh()
     79     mesh_B = V_B.mesh()

NotImplementedError: 
Sorry, only implemented for affine maps for now. To do non-affine, we'd need to
import much more of the assembly engine of UFL/TSFC/etc to do the assembly on
each supermesh cell.