Two Different PDE solvers give different answers

dajopez · December 3, 2020, 4:04am

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

load "SLEPc"
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.

prj · December 3, 2020, 6:43am

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?

dajopez · December 3, 2020, 12:32pm

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.

prj · December 3, 2020, 12:41pm

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.

dajopez · December 3, 2020, 1:10pm

Ok I’ll try that.

What does doing that do?

prj · December 3, 2020, 1:12pm

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

dajopez · December 3, 2020, 1:52pm

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

prj · December 3, 2020, 2:13pm

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.

dajopez · December 3, 2020, 2:40pm

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

dajopez · December 3, 2020, 2:41pm

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

load "SLEPc"
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;

dajopez · December 3, 2020, 2:41pm

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;

prj · December 3, 2020, 3:44pm

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

dajopez · December 3, 2020, 4:00pm

Like this?

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

prj · December 3, 2020, 4:01pm

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

dajopez · December 3, 2020, 4:04pm

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

prj · December 3, 2020, 4:06pm

Here is the link to the corresponding SLEPc documentation.

dajopez · January 12, 2021, 6:31am

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

prj · January 12, 2021, 6:36am

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

Topic		Replies	Views
Negative values in solution General Discussion	6	175	August 7, 2023
The elasticity problem with "varf" General Discussion	23	589	October 19, 2022
Doubt in the variational formulation of LaplaceRT.edp General Discussion	5	206	May 24, 2023
Cannot get example to work General Discussion	1	306	September 8, 2020
Problem with three equations General Discussion	1	325	December 12, 2020

Two Different PDE solvers give different answers

Related Topics