I have a problem while defining ratio of two solution functions. I am not getting the expected answer as seen in the pictures below. I would like to know whether I am understanding the issue properly or not. I fear, I won’t be able to post the exact code because it is too long but I am posting below a MWE which to my understanding is working perfectly fine(because there aren’t any zero values involved). I am also posting pictures below which are calculated in DG-0 space, to my understanding I expect zero values to the left and right of the non-zero layer of values(I am adding a tolerance value to actual denominator array to avoid issue of dividing by zero) , but that’s not what I get. I want to know whether the division is happening cell-wise for both solutions or not.
MWE
#calculating f expression
msh = create_unit_square(MPI.COMM_WORLD,1,1)
x = ufl.SpatialCoordinate(msh)
f = x[0] + x[1]
V = dolfinx.fem.functionspace(msh, ("DG", 0))
f_vl = fem.Function(V)
f_exp = dolfinx.fem.Expression(f,V.element.interpolation_points())
f_vl.interpolate(f_exp)
# calculating g expression
x = ufl.SpatialCoordinate(msh)
g = 2*(x[0])
V = dolfinx.fem.functionspace(msh, ("DG", 0))
g_vl = fem.Function(V)
g_exp = dolfinx.fem.Expression(g,V.element.interpolation_points())
g_vl.interpolate(g_exp)
# calculate ratio
ratio = fem.Function(V)
ratio.x.array[:] = f_vl.x.array/g_vl.x.array
dr @francesco-ballarin I appreciate your response, yes, these pictures don’t correspond to the code, but I am using similar approach to compute these. The issue is first solution is result of a very long code for which MWE is not reproducible, but that part is correct. I know I am calculating pic 1 and pic 2 correclty, my issue is only with ratio part. Please allow me to put the code how pic 2 is computed and for the ratio part.
# calculate pic2 i.e f- a.grad(u_h)
# compute gradient and interpolate into DG-0 space
mesh_coarse = u_H1_projection.function_space.mesh
gdim = mesh_coarse.geometry.dim
dg_space = dolfinx.fem.functionspace(mesh_coarse,('DG',0,(gdim,)))
grad_u_H1 = dolfinx.fem.Function(dg_space)
grad_uh = dolfinx.fem.Expression(grad(u_H1_projection),dg_space.element.interpolation_points())
grad_u_H1.interpolate(grad_uh)
# compute c_T as interpolation in DG-0 space
vec_f = Constant(mesh_coarse,ScalarType((np.cos(np.pi/3),np.sin(np.pi/3))))
f = Constant(mesh_coarse, ScalarType(0))
c_space = dolfinx.fem.functionspace(mesh_coarse,('DG',0))
c_T = dolfinx.fem.Function(c_space)
c_T_expr = dolfinx.fem.Expression(f - inner(vec_f,grad_u_H1),c_space.element.interpolation_points())
c_T.interpolate(c_T_expr)
# compute the ratio a_T/c_T approx tau
tol = 1e-8
c_T_new = dolfinx.fem.Function(c_space)
c_T_new.x.array[:] = c_T.x.array + tol
a_T_c_T = dolfinx.fem.Function(c_space)
a_T_c_T.x.array[:] = u_prime_avg.x.array/c_T_new.x.array
here in the ratio part u_prime_avg is calculated correctly which corresponds to pic 1.
However, let me note that adding 1e-8 isn’t a foolproof way of guaranteeing that the denominator is non-zero. Your pic2 shows that the solution can have negative values, so if by chance the Function you put at the denominator had value -1e-8 at one DOF you would still end up with a zero denominator.
dr @francesco-ballarin thank you for your feedback, my only issue is the ratio picture i,e the third pic, I want zero values on the left and right hand side of the layer(non-zero value layer) . I wonder whether code is doing corresponding cell-wise division of pic1 and pic2 or not. As, I can see non-zero values in the left bottom corner of pic3, seems division is not happening corresponding cell-wise( I might be wrong). Can you please provide insight on how can I guarantee here my denominator is not zero. I hope that might resolve my issue.
I’m not sure anything is actually wrong. You are potentially dividing two small values by each-other, and thus ratios can blow up (imagine two entries being of magnitude 1e-15, but one is 9e-15, the other one 9e-16).
