Real[int] from varf in PETSc

General Discussion
1

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)

2

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
3

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.

4

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

5

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?

6

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

7

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