This code produces usually the following plot:
However if i run the code multiple times, this plot sometimes appears:
if i wanted to use transfer of functions from one space to another iterative over time and have to run this code lets say 1000 times, this doesn’t seem very reliable.
Is there anything i can do to get the ‘same’ results all the time?
I obtain the same figure every time I run the program. You can also inspect the norm, as Im printing: print(func2dz.vector().norm("l2")) which to me is constant at: 16.392286467209278
Hey, i changed the code a little to make it visible how often the results vary, i used the center point of the result to compare the results. Try this code:
from fenics import *
import matplotlib.pyplot as plt
import math as mth
parameters['allow_extrapolation'] = True
height = 1.0
width = 1.0
length = 1.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length,width,height), 9, 9, 9)
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for f in facets(mesh):
marker[f] = height - DOLFIN_EPS < f.midpoint().z() < height + DOLFIN_EPS
submesh = MeshView.create(marker, 1)
V2d = FunctionSpace(submesh, 'CG', 2)
V3d = FunctionSpace(mesh, 'CG', 2)
func2d = project(Expression('x[0]*x[0] + x[1] * x[1]', degree = 5), V2d)
func3d = interpolate(func2d, V3d)
func3dz = func3d.dx(2)
n = 100
a = 0
b = 0
c = 0
for i in range (n):
func2dz = project(func3dz, V2d, solver_type = 'mumps')
centervalue = func2dz(0.5, 0.5, height)
# print(centervalue)
if centervalue == -0.6522767164072237:
a += 1
if centervalue != -0.6522767164072237:
b += 1
if mth.isnan(centervalue) == True:
c += 1
print(a)
print(b)
print(c)
plot(func2dz)
The value for the comparison was just the value i observed the most. However the The results are far from consistent, the last run, resulted in 92 times the expected value, 8 times a different value, from which 5 were nan.
I do not get the same result as you. Running your script with the docker container above, I got:
fenics@2c50fb5e1fbe:~/shared$ python3 test_projection.py
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Building point search tree to accelerate distance queries.
Computed bounding box tree with 323 nodes for 162 points.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
0
100
0
What is the value you get more regularly than the others and is this value always the same?
I did the same test with the interpolations, there i dont get different results, instead its always the same.
I see some of the variability you get as well. However, there are a few things I would like you to adress:
You are projecting the derivative of a function into the same space as your original function. This will not be a good representation, as the deriviatve of a CG-2 function is in a DG-1 space. See for instance: Interpolate the gradient of a dolfin function - #3 by VKTR
You should check that the interpolation of func2d into func3d gives a consistent result before trying to project the derivative.
Thats interesting, if i don’t compute the derivative i get the same results everytime.
Now how do i work with derivative? I could indeed use the suggested space, however another functionspace will not line up with the rest of the programm from which i extracted this problem.
Would a finite difference in the z-direction be adequate to workaround the problem?
I am not sure about your particular problem, but you should think about if it actually makes sense to project it into a CG 2 space when the function is DG-1. This may cause spurious behavior at a discontinuous boundary.
Try for instance projecting a function that is one in part of your domain and zero otherwise into CG1
from dolfin import *
mesh = UnitSquareMesh(10,10)
V = FunctionSpace(mesh, "CG", 2)
expr = Expression("x[0] < 0.5 ? x[0]-x[1] : 0", degree=3)
u = project(expr, V)
File("u.pvd") << u
File("u_p1.pvd") << project(Expression("x[0]-x[1]", degree=3), V)
Well the problem is a coupled equation of biharmonic equation on a plate, with heat equation in a volume. I solve the equation iterative and transfer the solutions by projecting or interpolating. However i am using the two functionspaces in which the equations and functionspaces are defined. The equations are also solved over time and are therefore solved e.g. a 1000 times. Now i am a bit unsure how a third space (DG-1) comes into place and how i would make the interfaces between the spaces.
The iterative approach was chosen by my supervisor in the very beginning of our work. I proposed a stronger coupled approach in a mixed space without any projection or interpolation. However we still wanted to try everything before chaning the approach. The current state is that interpolating and projecting actually works fine, but not with the derivative, what methods of deriving could i try? Are there any other ways to get the derivative in z direction?
I would suggest making your own projection or interpolation method. As the gradient of yourfunction is not uniquely defined at the degrees of freedom, you need to chose its given behavior.
Is it really an issue of parts of the domain? I am just a bit confused, since i just want the derivative in a certain direction, that is numerically stable for a number of repeated calculations.
I am just giving you the tools to do what you would like to do. The code in the «parts of domain» approach is a working example of a custom projection.
In your case, you are trying to project a derivative of a 3D field into 2D. There are several things you could consider:
-first project func3dz into an appropriate 3D space (DG degree-1)
then i think you should be able to interpolate into 2D.
You should Also inspect the behavior of func3dz after projection using xdmf.write_checkpoint and visualize it with paraview.
