I am currently working on solving the heat diffusion-advection equation on a 3D mesh created with gsmh, which includes several subdomains. In each subdomain of the mesh, I have to manually define a different fluid flow.
Initially, I defined a discontinuous cell-wise constant function to assign a fluid flow vector to each cell of the mesh (DG0), in a similar way to that used to define multiple materials:
Q = VectorFunctionSpace(mesh, ("DG", 0))
q = Function(Q)
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
block_size = Q.dofmap.index_map_bs
for i in range(num_cells):
# in subdomains 103, 104 and 105 the flow is null
if (cell_tags.values[i] == 103) | (cell_tags.values[i] == 104) | (cell_tags.values[i] == 105) :
q.x.array[[i*block_size, i*block_size+1, i*block_size+2]] = [0., 0., 0.]
# in subdomain 106 the flow is uniform
elif cell_tags.values[i] == 106:
q.x.array[[i*block_size, i*block_size+1, i*block_size+2]] = [0., 0., -1]
# in subdomain 107 the flow is radial
elif cell_tags.values[i] == 107:
# obtain coordinates of the centroid of the element
coord_centroid = Q.tabulate_dof_coordinates()[i]
# calculate unit vector in the direction of flux
u = coord_centroid[:2]/np.linalg.norm(coord_centroid[:2])
# assign the flow
q.x.array[[i*block_size, i*block_size+1, i*block_size+2]] = list(u) + [0.]
Now I would like to define the same fluid flow on continuous Galerkin elements of degree 2 (CG2). Do you have any ideas on how I could define the fluid flow on this new VectorFunctionSpace? Is this approach even possible? I tried to loop over the cells and then retrieve the degrees of freedom to assign the flow, but I am stuck at this stage:
Q = VectorFunctionSpace(mesh, ("CG", 2))
q = Function(Q)
num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
block_size = Q.dofmap.index_map_bs
for i in range(num_cells):
# get the DOF of the nodes of the cell
cell_dofs = Q.dofmap.cell_dofs(cell_tags.values[i])
for dof in cell_dofs:
q.x.array[...] = ...
Using a continuous element (“Langrange”, 2) to describe something that is cell-wise continuous is not well defined (A value at a vertex or an edge can take in multiple values).
If you really want to do something like this, I would suggest that you
Put the velocity in the DG, 0 space
Project the velocity into your second order continuous space.
What is the motivation for moving the velocity into a continuous space?
Thank you very much for your response @dokken.
In fact, in my advection-diffusion simulation, I impose a heat flux on one of the boundaries and a zero heat flux condition on all the other boundaries (so only Neumann conditions). Then I compute the energy given to the system versus the energy stored in the entire domain at each time step, and I noticed that they are different (both increase linearly but their slopes are different). Note that in my simulation, the Péclet number is close to 1 and I do not observe any oscillations or numerical instabilities.
So I suppose there is a problem in my simulation and I am exploring several possibilities, including the first one which would be related to the interpolation degrees of the different FunctionSpace.
In the example “A system of advection-diffusion-reaction equations” from the FEniCS tutorial, I saw that for the velocity field, CG elements of order 2 are used while for the concentrations CG elements of order 1 are used, and I was wondering if in my case it would be necessary to have a higher degree of interpolation for the velocity field compared to temperature.
That’s why I tried to have a velocity field on CG2 elements instead of DG0 elements.
Do you think this is not necessary for my simulation or that the problem may come from elsewhere?
Thanks again for all of your help.
The tutorial uses the velocity field from a previous tutorial, Where a P2-P1 element pair was used to solve for velocity and pressure (called Taylor-Hood), which is inf-sup stable.
Thank you for this clarification.
Does this mean that I can keep the flow defined on DG0 elements without it causing any issues for the diffusion-advection simulation?