I will probably post a demo at the scifem webpage of this soon, based of
2 Likes
Well, it took a while to figure out as I’d never used Docker before, but it turned out to be a lot simpler than I was making it out to be, I simply had to run the command line:
docker run --init -ti -p 8888:8888 dolfinx/lab:nightly
Now, before diving into it, I wanted to take a second to appreciate the paper you dug up, as it is actually quite heavily ingrained in a case study I have been working on, but have had little progress in.
With that out of the way, I needed to spend a few days mulling over the first few sections of the paper as these concepts are all still certainly foreign to me, so I apologize for taking a while to get back to you. It looks like you recreated the first example covered in section 5, and this is described under 5.1.1 titled “Div-conforming spaces”. There are still some questions I have:
- First, I wanted to confirm something about
It looks like you, in the code you provided initially, enforced tangency of \boldsymbol\sigma using the Lagrange multipliers following 5.1.2, and with this method, we instead solve for a system of now 3 equations:
\boldsymbol\sigma - \nabla u - l\textbf{k} = 0
\nabla \cdot \boldsymbol\sigma + r = g
\boldsymbol \sigma \cdot \textbf{k} = 0
where up is the zenith unit vector \textbf{k}, and w is the corresponding test function of \boldsymbol \sigma which in the paper was originally designated as \boldsymbol\tau, and the combination of these equations yields the equation that we are now solving for
<\boldsymbol\sigma, \boldsymbol\tau> - <\boldsymbol\tau,\nabla u> + <\boldsymbol\sigma,\nabla v> - <l,\boldsymbol\tau \cdot \textbf{k}> + <\gamma,\boldsymbol \sigma \cdot \textbf{k}> - <r,v> + <t,u> = -<g,v>.
-
I notice in the github post where you attempted to solve using H-div conforming spaces which reduces the problem to just solving for (38) and (39) via (44), that the L2 error for both
u_handsigma_hwere much higher than it was for the Lagrange multiplier method. Do you know if this due to limitations to how modern Fenics handles the Piola transform, or is this an advantage of using the Lagrange multiplier unrelated to Fenicsx? In the paper, it says that if we use a div-conforming space, vector fields are naturally mapped from a reference cell to each physical cell via the contravariant piola transform, so for this case, the sigma vectors should be more accurate if we expect Fenics to accurately handle this automatically. If you’re not certain, I feel it is worth the developers of Fenics to take a look at. -
To further test my second question, I was trying to modernize the
linear_shallow_water.py, and compare results to those inenergies.py(due to a lack of an exact solution as a function of time provided to compare with), but I find that I’m not well versed enough with Dolfinx to intuitively convert from legacy Dolfin to modern Dolfinx as I got stuck at rewritinginitial_conditions. Would you be able to test this case as well?
Oh and are you able to plot sigma? Should be a vector field, so using glyphs should work and we’d expect them to be tangential to the surface. I’m just having trouble on my end adding the sigma values to the grid.
Could you be a bit more specific here. What L^2 errors in my github post (and which github post), and what meshes have been used, as well as the order of the geometry. Are the examples you comparing using the same mesh and equivalent orders of discretization?
For instance, in: Add manifold example with real spaces · Issue #118 · scientificcomputing/scifem · GitHub
when using a third order geometry, I get small L2 errors:
L2 error in u: 6.449031074513135e-07
L2 error in sigma: 1.8679131149781764e-06
H(div) error in sigma: 1.408483169816556e-05
Here I do not use the original meshes by Marie, as they are only linear (affine) triangles, which can reduce the order of convergence, as discussed in her aforementioned paper.
We would usually expect these convergence rates to in-
crease by one order when we change the spaces to DG3
1 −
CG2. However, as discussed in (Bernard et al., 2008), higher-
order convergence can only be achieved if higher-order ap-
proximations to the manifold itself are used, and in our im-
plementation we use affine triangles. Hence, for DG3
1 − CG2,
we also observe second and first order convergence for L2
and H1 norms, respectively. This example tests the use of
three-dimensional vector fields on a two-dimensional mani-
fold mesh.
I might convert the linear shallow water equations eventually, but it is not high on my list of priorities.
THe initial conditions seems to be quite straightforward, as one is assigning zeros to all dofs, the other is a projection, which is for instance covered in: Where to find 'project'-function in dolfinx?
You would have to project sigma into a vector DG space for visualization (as VTK/Pyvista/Paraview only support arbitrary order Lagrange elements). Interpolation is currently broken from div-conforming spaces on manifolds, due to Finite element spaces on manifolds do not carry information about their value_shape in physical space · Issue #3619 · FEniCS/dolfinx · GitHub
1 Like
I’m sorry, my code was wrong. I was using the code from the bug report. What you have here is indeed highly accurate. Sorry for the misunderstanding! With error values for all solution variables being as low as they are, there is much less room for doubt.
That’s fine, I’ll continue tinkering with it in my own time. A lot of this stuff I find highly valuable for what I do, so the incentive is one-sided. Again, this was a great find!
The line:
u0 = interpolate(Constant((0.0,0.0,0.0)), S)
seems to mean that 0.0 is being assigned to 3 different dimensions, but the Dolfinx equivalent seems to only allow for 2 which is the same as the topological dimension of the mesh. The following compiled:
def expr(x):
init_values = np.zeros(x[:S.mesh.topology.dim].shape)
print(init_values.shape)
return init_values
u0 = fem.Function(S)
u0.interpolate(expr)
Hopefully that’s correct.
I’m just now looking at this, but in order to project the UFL expression onto the functionspace, it looks like I’m solving some variational problem, is that right? From the code, it looks like the variational form is
\int_\Omega f_h v dx = \int_\Omega fv dx
where v is the trial function and f is the UFL expression that we want to project onto V. I’ll have to give that a go, and hopefully it all works out.
I think I managed to get it to work without the interpolation step blowing up.
W = dolfinx.fem.functionspace(mesh, ("DG", 0, (mesh.geometry.dim, )))
sigma_d = dolfinx.fem.Function(W)
expr_d = dolfinx.fem.Expression(ufl.as_vector(sigma_h),W.element.interpolation_points)
sigma_d.interpolate(expr_d)
num_dofs = W.dofmap.index_map.size_local
num_verts = W.mesh.topology.index_map(0).size_local
values_d = np.zeros((num_dofs, 3), dtype=np.float64)
values_d[:, :mesh.geometry.dim] = sigma_d.x.array.real.reshape(num_dofs, W.dofmap.index_map_bs)
grid = pv.UnstructuredGrid(*dolfinx.plot.vtk_mesh(Ve))
grid["sigma"] = values_d
glyphs = grid.glyph(orient = "sigma", scale = "sigma", factor = 4)# tolerance=0.1)
grid.set_active_vectors("sigma")
u_plotter = pv.Plotter()
u_plotter.add_mesh(grid,show_edges= False, color = 'white')
u_plotter.add_mesh(glyphs, cmap = 'bwr')#, clim = [np.min(grid["sigma"]),np.max(grid["sigma"])])
u_plotter.view_xy()
#u_plotter.show()
file = u_plotter.screenshot("my_pyvista_screenshot.png")
and compare that to the sigma plot from the paper:
so it looks like it works. I asked because I was having some trouble last night, and I also had seen your issue beforehand thinking that might have something to do with it, but I think I have it working now.
1 Like