 # Two Different PDE solvers give different answers

Hello, I am trying to switch my code to parallel, but the new solver I use gives out wrong data. I know its wrong because it doesnt line up with analytical results. Here is the correct implementation:

``````// Problem
// OP = A - sigma B ; // the shifted matrix
varf op (u1, u2)
= int2d(Th)(
dx(u1)*dx(u2)
+ dy(u1)*dy(u2)
- sigma* u1*u2
)
+ on(a0,a1,a2,a3,a4,a5, u1=0)
;

varf b ([u1], [u2]) = int2d(Th)(u1*u2); //no boundary condition

matrix OP = op(Vh, Vh, solver=Crout, factorize=1); //crout solver because the matrix in not positive
matrix B = b(Vh, Vh, solver=CG, eps=1e-20);

int k = EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV,
tol=1e-10, maxit=0, ncv=0);

``````

Here is the code I cant get to work:

``````
// Problem
// OP = A - sigma B ; // the shifted matrix
varf op (u1, u2)
= int2d(Th)(
dx(u1)*dx(u2)
+ dy(u1)*dy(u2)
- sigma* u1*u2
)
+ on(a0,a1,a2,a3,a4,a5, u1=0)
;

varf b ([u1], [u2]) = int2d(Th)(u1*u2); //no boundary condition

Mat OP;
include "macro_ddm.idp"
createMat(Th, OP, Pk)
OP = op(Vh, Vh, tgv = -2,sym=1); //crout solver because the matrix in not positive
matrix loc = b(Vh, Vh, tgv = -20, sym =1);
Mat B(OP, loc);

int k = EPSSolve(OP, B, values=ev, vectors=eV,sparams = "-st_type sinvert -eps_nev " + nev + " -eps_target " + sigma);

``````

Any direction will be greatly appreciated.

As I’ve already told you, put your code between three single quote ```, otherwise it is not copy/paste-able. Make sure it’s runnable, it’s clearly not the case here, e.g., what’s `sigma`?

``````code should appear like this;
``````

It’s wrong to put the shift in the left-hand side matrix. You have to remove `- sigma u1 * u2` from the varf.

How is this post different than this one? Also, this question is not related to PDE solvers, but eigensolvers, if I’m not mistaken? What are the eigenvalues you are looking for? Closest to the shift? What are the “analytical results” you are mentioning, i.e., exact values you are supposed to get?

My single quotes are different than yours, they are slanted the other way and it won’t take. I had to copy and paste from mobile.

So I want to solve Laplace(u) = lambda u with u= 0 on the boundary. I know specific things about solutions to this on the boundary that are independent of eigenvalue or the shape of the triangle. Solutions align with this results in the case that’s working vs the case that’s not. These values are also computed independent of the the solution so different solvers won’t affect them at all. I can’t see why there is a discrepancy.

Yeah I just got confused I want to use the eigensolver. Sigma is always 1 in my code. The solvers also put out the same eogenvalues but give out different functions.

Again, you need to remove `- sigma* u1*u2` from your `op` varf.
You are also missing `+ on(a0,a1,a2,a3,a4,a5, u1=0)` in your `b` varf.

Ok I’ll try that.

What does doing that do?

1. no shift in the varf, because you shouldn’t have to deal with that yourself, it’s handled by SLEPc internally
2. use the correct boundary conditions

I made those changes and it gives the exact same results as before

I can’t run either of your script, they are not self-contained: what is `Th`, `Vh`… ?
If you put copy/paste-able runnable scripts, I can have a look.

Let me post the full scripts so you can look. I cut out segments because a lot of it is just making arrays and calculating values

Parallel one that doesnt work

``````// Parameters
real sigma = 1; //value of the shift (ignore)
real a = 1./sqrt(3);   //length of bottom side
string lengthName = "306090";
int n = 100;   //numer of mesh points on each side
int nev = 100; //number of computed eigen values

//Important calculated values
real triArea = .5*a; //area of triangle

real leftLength = 1;
real bottomLength = a;
real hypotLength = sqrt(leftLength^2. + bottomLength^2.);

real leftNeumann = leftLength/triArea;      // Exact Neumann data value of left side
real bottomNeumann = bottomLength/triArea;  // Exact Neumann data value of bottom side
real hypotNeumann = hypotLength/triArea;    // Exact Neumann data value of hypotenuse

// Mesh
border a0(t=0., .5){x=0; y=1.-t;}        //top half of left side
border a1(t=0., .5){x=0; y=.5-t;}        //bottom half of left side
border a2(t=0., .5){x=a*(1.-t); y=t;}    //bottom half of hypotenuse
border a3(t=0., .5){x= a*(.5-t); y=.5+t;}//top half of hypotenuse
border a4(t=0., .5){x=a*t; y=0.;}        //left half of bottom side
border a5(t=0., .5){x=a*(.5+t); y=0.;}   //right half of bottom side

//Freefem stuff. Builds mesh and declares function space. Stuff in this section from example.
mesh Th = buildmesh(a0(n) + a1(n) + a2(n)+ a3(n)+ a4(n)+ a5(n));
func Pk = P2;
fespace Vh(Th, Pk);
Vh u1, u2;

// Problem
// OP = A - sigma B ; // the shifted matrix
varf op (u1, u2)
= int2d(Th)(
dx(u1)*dx(u2)
+ dy(u1)*dy(u2)
)
+ on(a0,a1,a2,a3,a4,a5, u1=0)
;

varf b ([u1], [u2]) = int2d(Th)(u1*u2)
+ on(a0,a1,a2,a3,a4,a5, u1=0); //no boundary condition

Mat OP;
include "macro_ddm.idp"
createMat(Th, OP, Pk)
OP = op(Vh, Vh, tgv = -2,sym=1); //crout solver because the matrix in not positive
matrix loc = b(Vh, Vh, tgv = -20, sym =1);
Mat B(OP, loc);

// important remark:
// the boundary condition is make with exact penalization:
// we put 1e30=tgv on the diagonal term of the lock degree of freedom.
// So take Dirichlet boundary condition just on \$a\$ variational form
// and not on \$b\$ variational form.
// because we solve \$ w=OP^-1*B*v \$

// Solve
real[int] ev(nev); //to store the nev eigenvalue

real[int] accuracyLeft(nev); //stores accuracy of left face (%)
real[int] accuracyBottom(nev); //stores accuracy of bottom face (%)
real[int] accuracyHypot(nev); //stores accuracy of hypotenuse (%)

real[int] ratioLeft(nev); //stores ratio of two halves of left face
real[int] ratioBottom(nev); //stores ratio of two halves of Bottom face
real[int] ratioHypot(nev); //stores ratio of two halves of hypotenuse

real[int] computedNeumannBottomLeftList(nev);
real[int] computedNeumannBottomRightList(nev);
real[int] computedNeumannLeftTopList(nev);
real[int] computedNeumannLeftBottomList(nev);
real[int] computedNeumannHypotTopList(nev);
real[int] computedNeumannHypotBottomList(nev);

real[int] volumeDxList(nev);
real[int] volumeDyList(nev);

Vh[int] eV(nev); //to store the nev eigenvector

int k = EPSSolve(OP, B, values=ev, vectors=eV,sparams = "-st_type sinvert -eps_nev " + nev + " -eps_target " + sigma);

// Display & Plot
for (int i = 0; i < nev; i++){
u1 = eV[i];

//Calculated Neuman Boundary data
real computedNeumannBottomLeft = int1d(Th,a4)(dy(u1)^2)/ev[i];
real computedNeumannBottomRight = int1d(Th,a5)(dy(u1)^2)/ev[i];

real computedNeumannLeftTop = int1d(Th,a0)(dx(u1)^2)/ev[i];
real computedNeumannLeftBottom = int1d(Th,a1)(dx(u1)^2)/ev[i];

real computedNeumannHypotTop = (1./(a^2.+1.))*int1d(Th,a3)((a*dy(u1)+dx(u1))^2)/ev[i];
real computedNeumannHypotBottom = (1./(a^2.+1.))*int1d(Th,a2)((a*dy(u1)+dx(u1))^2)/ev[i];

real volumeDx = int2d(Th)(dx(u1)*dx(u1))/ev[i];
real volumeDy = int2d(Th)((dy(u1)^2))/ev[i];

//Stores Neumann Data
computedNeumannBottomLeftList[i] = computedNeumannBottomLeft;
computedNeumannBottomRightList[i] = computedNeumannBottomRight;
computedNeumannLeftTopList[i] = computedNeumannLeftTop;
computedNeumannLeftBottomList[i] = computedNeumannLeftBottom;
computedNeumannHypotTopList[i] = computedNeumannHypotTop;
computedNeumannHypotBottomList[i] = computedNeumannHypotBottom;
volumeDxList[i] = volumeDx;
volumeDyList[i] = volumeDy;

//Stores percentage relative accuracy of Neumann Data all various faces
accuracyLeft[i] = (computedNeumannLeftTop + computedNeumannLeftBottom - leftNeumann)*100./leftNeumann;
accuracyBottom[i] = (computedNeumannBottomLeft+computedNeumannBottomRight - bottomNeumann)*100./bottomNeumann;
accuracyHypot[i] = (computedNeumannHypotTop+computedNeumannHypotBottom-hypotNeumann)*100./hypotNeumann;

//Stores ratio of two halves of the same face. One touching hypotenuse is the numerator. For hypot, the

ratioLeft[i] = computedNeumannLeftTop/computedNeumannLeftBottom;
ratioBottom[i] = computedNeumannBottomRight/computedNeumannBottomLeft;
ratioHypot[i] = computedNeumannHypotTop/computedNeumannHypotBottom;

//real gg = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1));
//real mm = int2d(Th)(u1*u1) ;
//cout << "lambda[" << i << "] = " << ev[i] << ", err= " << int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) - (ev[i])*u1*u1) << endl;
//plot(eV[i], cmm="Eigen Vector "+i+" value ="+ev[i], wait=true, value=true);
}

cout << "L2mass: "<< int2d(Th)(u1^2)^.5 << endl;
cout << "Neumann Data Bottom Left Face: " << computedNeumannBottomLeftList << endl;
cout << "Neumann Data Left Bottom Face: " << computedNeumannLeftBottomList << endl;
cout << "Neumann Data Bottom Right Face: " << computedNeumannBottomRightList << endl;

cout << "Neumann Data Left Top Face: " << computedNeumannLeftTopList << endl;
cout << "Neumann Data Diagonal Bottom Face: " << computedNeumannHypotBottomList<< endl;
cout << "Neumann Data Diagonal Top Face: " << computedNeumannHypotTopList << endl;

cout << "Accuracy for Left face (%): " << accuracyLeft << endl;
cout << "Accuracy for Bottom face (%): " << accuracyBottom << endl;
cout << "Accuracy for Hypotenuse (%): " << accuracyHypot << endl;

cout << "Ratio for Left face: " << ratioLeft << endl;
cout << "Ratio for Bottom face: " << ratioBottom << endl;
cout << "Ratio for Hypotenuse: " << ratioHypot << endl;

cout << "Volume dx integral: " << volumeDxList << endl;
cout << "Volume dy integral: " << volumeDyList << endl;

//File Saving
ofstream file("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannBottomLeft.txt");
file << computedNeumannBottomLeftList << endl;

ofstream file2("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannBottomRight.txt");
file2 << computedNeumannBottomRightList << endl;

ofstream file3("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannLeftBottom.txt");
file3 << computedNeumannLeftBottomList << endl;

ofstream file4("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannLeftTop.txt");
file4 << computedNeumannLeftTopList << endl;

ofstream file5("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannHypotBottom.txt");
file5 << computedNeumannHypotBottomList << endl;

ofstream file6("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannHypotTop.txt");
file6 << computedNeumannHypotTopList << endl;

ofstream file7("length" + lengthName + "nodes" + n + "eVals" + nev + "accuracyLeft.txt");
file7 << accuracyLeft << endl;

ofstream file8("length" +lengthName + "nodes" + n + "eVals" + nev + "accuracyBottom.txt");
file8 << accuracyBottom << endl;

ofstream file9("length" + lengthName + "nodes" + n + "eVals" + nev + "accuracyHypot.txt");
file9 << accuracyHypot << endl;

ofstream file10("length" + lengthName + "nodes" + n + "eVals" + nev + "ratioLeft.txt");
file10 << ratioLeft << endl;

ofstream file11("length" + lengthName + "nodes" + n + "eVals" + nev + "ratioBottom.txt");
file11 << ratioBottom << endl;

ofstream file12("length" + lengthName + "nodes" + n + "eVals" + nev + "ratioHypot.txt");
file12 << ratioHypot << endl;

ofstream file13("length" + lengthName + "nodes" + n + "eVals" + nev + "VolumeIntegralDx.txt");
file13 << volumeDxList << endl;

ofstream file14("length" + lengthName + "nodes" + n + "eVals" + nev + "VolumeIntegralDy.txt");
file14 << volumeDyList << endl;
``````

One I wrote that works

``````// Parameters
real sigma = 1; //value of the shift (ignore)
real a = 1./sqrt(3);   //length of bottom side
string lengthName = "306090";
int n = 350;   //numer of mesh points on each side
int nev = 1000; //number of computed eigen values

//Important calculated values
real triArea = .5*a; //area of triangle

real leftLength = 1;
real bottomLength = a;
real hypotLength = sqrt(leftLength^2. + bottomLength^2.);

real leftNeumann = leftLength/triArea;      // Exact Neumann data value of left side
real bottomNeumann = bottomLength/triArea;  // Exact Neumann data value of bottom side
real hypotNeumann = hypotLength/triArea;    // Exact Neumann data value of hypotenuse

// Mesh
border a0(t=0., .5){x=0; y=1.-t;}        //top half of left side
border a1(t=0., .5){x=0; y=.5-t;}        //bottom half of left side
border a2(t=0., .5){x=a*(1.-t); y=t;}    //bottom half of hypotenuse
border a3(t=0., .5){x= a*(.5-t); y=.5+t;}//top half of hypotenuse
border a4(t=0., .5){x=a*t; y=0.;}        //left half of bottom side
border a5(t=0., .5){x=a*(.5+t); y=0.;}   //right half of bottom side

//Freefem stuff. Builds mesh and declares function space. Stuff in this section from example.
mesh Th = buildmesh(a0(n) + a1(n) + a2(n)+ a3(n)+ a4(n)+ a5(n));
fespace Vh(Th, P2);
Vh u1, u2;

// Problem
// OP = A - sigma B ; // the shifted matrix
varf op (u1, u2)
= int2d(Th)(
dx(u1)*dx(u2)
+ dy(u1)*dy(u2)
- sigma* u1*u2
)
+ on(a0,a1,a2,a3,a4,a5, u1=0)
;

varf b ([u1], [u2]) = int2d(Th)(u1*u2); //no boundary condition

matrix OP = op(Vh, Vh, solver=Crout, factorize=1); //crout solver because the matrix in not positive
matrix B = b(Vh, Vh, solver=CG, eps=1e-20);

// important remark:
// the boundary condition is make with exact penalization:
// we put 1e30=tgv on the diagonal term of the lock degree of freedom.
// So take Dirichlet boundary condition just on \$a\$ variational form
// and not on \$b\$ variational form.
// because we solve \$ w=OP^-1*B*v \$

// Solve
real[int] ev(nev); //to store the nev eigenvalue

real[int] accuracyLeft(nev); //stores accuracy of left face (%)
real[int] accuracyBottom(nev); //stores accuracy of bottom face (%)
real[int] accuracyHypot(nev); //stores accuracy of hypotenuse (%)

real[int] ratioLeft(nev); //stores ratio of two halves of left face
real[int] ratioBottom(nev); //stores ratio of two halves of Bottom face
real[int] ratioHypot(nev); //stores ratio of two halves of hypotenuse

real[int] computedNeumannBottomLeftList(nev);
real[int] computedNeumannBottomRightList(nev);
real[int] computedNeumannLeftTopList(nev);
real[int] computedNeumannLeftBottomList(nev);
real[int] computedNeumannHypotTopList(nev);
real[int] computedNeumannHypotBottomList(nev);

real[int] volumeDxList(nev);
real[int] volumeDyList(nev);

Vh[int] eV(nev); //to store the nev eigenvector

int k = EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV,
tol=1e-10, maxit=0, ncv=0);

// Display & Plot
for (int i = 0; i < k; i++){
u1 = eV[i];

//Calculated Neuman Boundary data
real computedNeumannBottomLeft = int1d(Th,a4)(dy(u1)^2)/ev[i];
real computedNeumannBottomRight = int1d(Th,a5)(dy(u1)^2)/ev[i];

real computedNeumannLeftTop = int1d(Th,a0)(dx(u1)^2)/ev[i];
real computedNeumannLeftBottom = int1d(Th,a1)(dx(u1)^2)/ev[i];

real computedNeumannHypotTop = (1./(a^2.+1.))*int1d(Th,a3)((a*dy(u1)+dx(u1))^2)/ev[i];
real computedNeumannHypotBottom = (1./(a^2.+1.))*int1d(Th,a2)((a*dy(u1)+dx(u1))^2)/ev[i];

real volumeDx = int2d(Th)(dx(u1)*dx(u1))/ev[i];
real volumeDy = int2d(Th)((dy(u1)^2))/ev[i];

//Stores Neumann Data
computedNeumannBottomLeftList[i] = computedNeumannBottomLeft;
computedNeumannBottomRightList[i] = computedNeumannBottomRight;
computedNeumannLeftTopList[i] = computedNeumannLeftTop;
computedNeumannLeftBottomList[i] = computedNeumannLeftBottom;
computedNeumannHypotTopList[i] = computedNeumannHypotTop;
computedNeumannHypotBottomList[i] = computedNeumannHypotBottom;
volumeDxList[i] = volumeDx;
volumeDyList[i] = volumeDy;

//Stores percentage relative accuracy of Neumann Data all various faces
accuracyLeft[i] = (computedNeumannLeftTop + computedNeumannLeftBottom - leftNeumann)*100./leftNeumann;
accuracyBottom[i] = (computedNeumannBottomLeft+computedNeumannBottomRight - bottomNeumann)*100./bottomNeumann;
accuracyHypot[i] = (computedNeumannHypotTop+computedNeumannHypotBottom-hypotNeumann)*100./hypotNeumann;

//Stores ratio of two halves of the same face. One touching hypotenuse is the numerator. For hypot, the

ratioLeft[i] = computedNeumannLeftTop/computedNeumannLeftBottom;
ratioBottom[i] = computedNeumannBottomRight/computedNeumannBottomLeft;
ratioHypot[i] = computedNeumannHypotTop/computedNeumannHypotBottom;

//real gg = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1));
//real mm = int2d(Th)(u1*u1) ;
//cout << "lambda[" << i << "] = " << ev[i] << ", err= " << int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) - (ev[i])*u1*u1) << endl;
//plot(eV[i], cmm="Eigen Vector "+i+" value ="+ev[i], wait=true, value=true);
}

cout << "L2mass: "<< int2d(Th)(u1^2)^.5 << endl;
cout << "Neumann Data Bottom Left Face: " << computedNeumannBottomLeftList << endl;
cout << "Neumann Data Left Bottom Face: " << computedNeumannLeftBottomList << endl;
cout << "Neumann Data Bottom Right Face: " << computedNeumannBottomRightList << endl;

cout << "Neumann Data Left Top Face: " << computedNeumannLeftTopList << endl;
cout << "Neumann Data Diagonal Bottom Face: " << computedNeumannHypotBottomList<< endl;
cout << "Neumann Data Diagonal Top Face: " << computedNeumannHypotTopList << endl;

cout << "Accuracy for Left face (%): " << accuracyLeft << endl;
cout << "Accuracy for Bottom face (%): " << accuracyBottom << endl;
cout << "Accuracy for Hypotenuse (%): " << accuracyHypot << endl;

cout << "Ratio for Left face: " << ratioLeft << endl;
cout << "Ratio for Bottom face: " << ratioBottom << endl;
cout << "Ratio for Hypotenuse: " << ratioHypot << endl;

cout << "Volume dx integral: " << volumeDxList << endl;
cout << "Volume dy integral: " << volumeDyList << endl;

//File Saving
ofstream file("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannBottomLeft.txt");
file << computedNeumannBottomLeftList << endl;

ofstream file2("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannBottomRight.txt");
file2 << computedNeumannBottomRightList << endl;

ofstream file3("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannLeftBottom.txt");
file3 << computedNeumannLeftBottomList << endl;

ofstream file4("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannLeftTop.txt");
file4 << computedNeumannLeftTopList << endl;

ofstream file5("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannHypotBottom.txt");
file5 << computedNeumannHypotBottomList << endl;

ofstream file6("length" + lengthName + "nodes" + n + "eVals" + nev + "computedNeumannHypotTop.txt");
file6 << computedNeumannHypotTopList << endl;

ofstream file7("length" + lengthName + "nodes" + n + "eVals" + nev + "accuracyLeft.txt");
file7 << accuracyLeft << endl;

ofstream file8("length" +lengthName + "nodes" + n + "eVals" + nev + "accuracyBottom.txt");
file8 << accuracyBottom << endl;

ofstream file9("length" + lengthName + "nodes" + n + "eVals" + nev + "accuracyHypot.txt");
file9 << accuracyHypot << endl;

ofstream file10("length" + lengthName + "nodes" + n + "eVals" + nev + "ratioLeft.txt");
file10 << ratioLeft << endl;

ofstream file11("length" + lengthName + "nodes" + n + "eVals" + nev + "ratioBottom.txt");
file11 << ratioBottom << endl;

ofstream file12("length" + lengthName + "nodes" + n + "eVals" + nev + "ratioHypot.txt");
file12 << ratioHypot << endl;

ofstream file13("length" + lengthName + "nodes" + n + "eVals" + nev + "VolumeIntegralDx.txt");
file13 << volumeDxList << endl;

ofstream file14("length" + lengthName + "nodes" + n + "eVals" + nev + "VolumeIntegralDy.txt");
file14 << volumeDyList << endl;

``````

Can you try to add `-eps_gen_hermitian` to your `EPSSolve` `sparams`, please?

Like this?

``````int k = EPSSolve(OP, B, values=ev, vectors=eV,sparams = "-st_type sinvert -eps_nev " + nev + " -eps_target " + sigma + "-eps_gen_hermitian");
``````

No… you are missing a space, this will be parsed as `-eps_target 1-eps_gen_hermitian`.

That seems to have fixed it. Could you explain what those parameters mean?

Here is the link to the corresponding SLEPc documentation.

Sorry to bother again, but where is the ‘include macro_ddm.idp’ is coming from? I haven’t been able to find it in the SLEPc documentation

It’s a FreeFEM file, you should not have to edit this file.