# Real[int] from varf in PETSc

Hi FreeFEM community,

I have noticed some strange behavior when I use a `varf` to create a `real[int]` in PETSc, and I’m not sure if this is a bug or my own misunderstanding. Please see that attached MWE (I’ve added some explainer text that will print via cout when the script is run). Note that the discrepancies referred to in this script are tiny, but in my actual application, the discrepancies are resulting in a division by 0. I imagine this must have something to do with ghost cell values? Any insight will be helpful!

Cheers,
Chris

test.edp (2.2 KB)

If I use `-pc_type lu`, I get:

``````\$ ff-mpirun -n 4 test.edp -v 0
CASE 1: serial (w/o PETSc)
value 0 = 0.41000000000000031
value 1 = 0.41000000000000031
value 2 = 0.41000000000000031
All values are equal.
CASE 2: with PETSc
value 0 = 0.41000000000007553
value 1 = 0.41000000000000031
value 2 = 0.41000000000000031
I understand why value 0 and value 1 may not be equal...
But why does value 1 not equal value 2 in case 2 (even with 1 proc)?
``````

So the values are equal up to machine epsilon.

1 Like

Thanks, Pierre. I will try to come up with a better MWE this afternoon. In my problem, I was using a direct solver.

In the meantime, can you see any reason why a ` real[int]` constructed from `varf1`(where the test function is written into the form AND BC) might not be equal to one built from `varf2` (where the test function only appears in the BC)? Of course, I’m referring to a case where the test function does not change between these two situations.

I’m not sure I understand your question, as both `varf` in your code return the same vector.

My apologies, let me clarify with a better MWE:
test2.edp (1.0 KB)

The main part of interest is here:

``````varf varf1( def(u2), def(v) ) = int2d(Th) ( f(u1,u2,v) ) + on(inds, u2 = 0);
varf varf2( def(u3), def(v) ) = int2d(Th) ( f(u1,u2,v) ) + on(inds, u3 = 0);
A = varf1(Vh, Vh, tgv = -10);
real[int] vec0 = A*u2[];
real[int] vec1 = varf1(0, Vh, tgv=-10);
real[int] vec2 = varf2(0, Vh, tgv=-10);
``````

I would expect all of these to `vec` to be equal (within machine precision). However, I find that, while `vec0` and `vec2` are equal, `vec1` is different (it is zero). Am I doing something wrong/is that expected?

Your `varf1` has no linear part, only bilinear. So if you assemble a vector out of it, you get nothing.

Indeed, this was a pretty fundamental misunderstanding on my part. Thank you for clarifying!