Hello everyone,
I am trying to calculate the convergence rate of my method. And there is no exact solution, I wanted to have the finer mesh and smaller time step as the “exact solution”, and I run the test on the coarser mesh, then I can calculate the convergence rate.
The difficulty I have is, how to interpolate the solution on the finer mesh to the coarser, (I know how to interpolate with analytic solution but not this one).
My dolfinx version is 0.6.0
Thank you very much!
Best,
Rui
In V0.6.0, there is an interpolation operator between meshes: Interpolation between different meshes by massimiliano-leoni · Pull Request #2245 · FEniCS/dolfinx · GitHub
See: https://github.com/FEniCS/dolfinx/blob/v0.6.0/python/test/unit/fem/test_function.py#L177-L218
Note that this operator has been improved in release v0.7.2, to work better with meshes with non-aligning outer boundaries.
how to interpolate the solution on the finer mesh to the coarser
Also note (perhaps it’s a typo) you should interpolate from the coarse mesh to the fine mesh for error measurement; as your coarse mesh cannot resolve the fine solution.
If your solutions are smooth, you could also employ a higher order basis function in your high resolution reference approximation. Then you don’t need to worry about coarse and fine meshes.
Thank you. When i try a simple 2d square mesh. I want to double the mesh point per side. and it gives me error. And sorry my version on the virtual environment is Dolfinx version: 0.5.2.
Traceback (most recent call last):
File “finer-mesh.py”, line 59, in
u1.interpolate(u0)
File “miniconda/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/function.py”, line 343, in interpolate
_interpolate(u, cells)
File “miniconda/envs/fenicsx-env/lib/python3.10/functools.py”, line 889, in wrapper
return dispatch(args[0].class)(*args, **kw)
File “miniconda/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/function.py”, line 324, in _
self._cpp_object.interpolate(u._cpp_object, cells)
RuntimeError: Interpolation on different meshes not supported (yet).
The below is the code:
#import cffi
import numpy as np
#import pytest
import ufl
from dolfinx.fem import (Function, FunctionSpace, TensorFunctionSpace,
VectorFunctionSpace, assemble_scalar, form)
from dolfinx.geometry import (BoundingBoxTree, compute_colliding_cells,
compute_collisions)
from dolfinx.mesh import (CellType, create_mesh, create_unit_cube,
create_unit_square, locate_entities_boundary,
meshtags)
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
from dolfinx import fem, mesh, io
from dolfinx.fem import Constant
from dolfinx.fem.assemble import assemble_matrix, create_vector, assemble_vector, apply_lifting, set_bc, assemble_scalar
import ufl
**n_nodes = 10 **
nx = n_nodes
ny = n_nodes
**mesh0 = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-1.0, -1.0]), np.array([1.0, 1.0])], [nx, ny], dolfinx.cpp.mesh.CellType.triangle) **
mesh1 = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([-1.0, -1.0]), np.array([1.0, 1.0])], [2nx, 2ny], dolfinx.cpp.mesh.CellType.triangle)
def f(x):
return (7 * x[1], x[2] + 0.4)
el0 = ufl.VectorElement(“Lagrange”, mesh0.ufl_cell(), 1, dim=2)
V0 = FunctionSpace(mesh0, el0)
el1 = ufl.VectorElement(“Lagrange”, mesh1.ufl_cell(), 1, dim=2)
V1 = FunctionSpace(mesh1, el1)
u0 = Function(V0)
u0.interpolate(f)
u0.x.scatter_forward()
u1 = Function(V1)
u1.interpolate(u0)
u1.x.scatter_forward()
u1_ex = Function(V1)
u1_ex.interpolate(f)
u1_ex.x.scatter_forward()
assert np.allclose(u1_ex.x.array, u1.x.array)
Interpolation between meshes was introduced in 0.6.0. To use this functionality, please update your environment.
Thank you, Nate! I have a question. If I use the Taylor-hood element P3-P2 vs P2-P1. When I calculate the rate
rate = log(error of (P2-P1)/ error(P2-P2))/ log(3/2) right? Since P3 is cubic, and on each edge, the leg is split into 3 parts, while in P2, the leg is split into 2 parts.
I’m not sure what you mean by the “leg” terminology.
Perhaps this script will help: DG Stokes demo.
I compute the numerical error in the L_2 norm where the approximation is is sought from a P2-P1 DG Taylor-Hood space whilst the “true solution” is approximated by interpolating onto the higher order P4-P3 DG Taylor-Hood space.
Thank you so much, Nate! My analytical solution is unknown.
I am hoping to have a mesh and calculate the rate by using different orders of Taylor hood space.
Let the solution on the higher-order Taylor-hood space say P5-P4 on the mesh be the “true solution”, and lower-order Taylor-hood P2-P1, P3-P2, and P4-P3 on the mesh be the coarse solution.
Take a 2D triangle element as an example. P2 means quadratic, so on each edge of the triangle, there will be 3 points including the endpoints. P3 means cubic, so on each edge of the triangle, there will be four points.
I need the ratio of the “h” to calculate the rate. I think the relationship of the “h” will be the number of segments of P3 over the number of segments of P2.
Thus:
rate = log(error of (P2-P1)/ error(P3-P2))/ log(3/2)
Am I understanding you correctly?
" If your solutions are smooth, you could also employ a higher order basis function in your high resolution reference approximation. Then you don’t need to worry about coarse and fine meshes."
Another thought I have: I can have the length of my solution vector, which is the degree of freedom.
I can use that ratio.
Thank you so much!
This depends on your definition of h. Typically it denotes a measure of the mesh cell size. A common measure is the cell diameter. This is useful in the context of hp-error analysis such that we may distinguish the error contribution from the spatial and spectral approximation. Note here that h is independent of p and vice versa. Therefore you shouldn’t be concerned with spatial distances between degrees of freedom. I worry that you’re mixing up h and p and will not get a reasonable estimate of the error you’re accumulating in your problem.