Poor results with PETSc

General Discussion
1

Hello,
I implemented a transient 3D piezoanalysis that works fine with problem. Now, I want to accelerate the simulation by using a parallel solver. However, I can’t get proper results using PETSc. The results seem to not converge properly, as I keep getting a very distorted mesh. I tried different solvers like lu or gmres. Also a “static” simulation (finding displacement U) does converge okay’isch.

Only the results in the first time step are accuarte (with lu). Could someone please help me with this problem?

load "msh3"
load "lapack"
load "PETSc"

//...

varf piezoTrans([phi, Fx,Fy,Fz], [Vphi, Vx, Vy, Vz]) =
	 int3d(Th)( dings(sigma(Ux,Uy,Uz), [Vx,Vy,Vz])+ dings2(D1(Ux,Uy,Uz),Vphi) )  
	 +int3d(Th)( dings2(D2(phi),Vphi) - dings(sigmaE(phi), [Vx,Vy,Vz]))
			  
	
	+int3d(Th)(Fx*Vx + Fy*Vy + Fz*Vz)
	
	 +on(elektrode1, phi=50)
	 +on(elektrode2, phi=0);

//this is exactly the same as above, but works
Vh phi2, Fx2,Fy2,Fz2, Vphi, Vx, Vy, Vz;
problem piezoTrans2([phi2, Fx2,Fy2,Fz2], [Vphi, Vx, Vy, Vz]) =
	 int3d(Th)( dings(sigma(Ux,Uy,Uz), [Vx,Vy,Vz])+ dings2(D1(Ux,Uy,Uz),Vphi) )  
	 +int3d(Th)( dings2(D2(phi2),Vphi) - dings(sigmaE(phi2), [Vx,Vy,Vz]))
			  
	
	+int3d(Th)(Fx2*Vx + Fy2*Vy + Fz2*Vz)
	
	 +on(elektrode1, phi2=50)
	 +on(elektrode2, phi2=0);

fespace Vh2(Th, [P1,P1,P1,P1]);

real dt = 20e-9;
real T = 10e-6;
int m, M = T/dt;

Vh Ux, Uy, Uz;

mesh3[int] anim(M);

matrix AA = piezoTrans(Vh2,Vh2, tgv=-1);
real[int] rhs = piezoTrans(0,Vh2, tgv=-1);
Mat A(AA);


for(m = 0; m < 3; m++){	
	AA = piezoTrans(Vh2,Vh2, tgv=-1);
	A = Mat(A);
	rhs = piezoTrans(0,Vh2, tgv=-1);
	
	
	
	set(A, sparams="-pc_type lu -ksp_view -ksp_monitor -ksp_monitor_true_residual -ksp_view_final_residual");
	//set(A, sparams="-ksp_view -ksp_monitor -ksp_monitor_true_residual -ksp_view_final_residual -pc_type gamg -ksp_type gmres");
	//set(A, sparams="-ksp_view -ksp_monitor -ksp_monitor_true_residual -ksp_view_final_residual -pc_type lu -ksp_type preonly -ksp_max_it 200 -n 4");
	
	Vh2 [phi,Fx,Fy,Fz];
	phi[] = A^-1*rhs;
	
	piezoTrans2;
	
	Vh deltaX, deltaY, deltaZ;
	
	deltaX = Fx-Fx2;
	deltaY = Fy-Fy2;
	deltaZ = Fz-Fz2;
	
	cout << "AbsX: " << (deltaX[]'*deltaX[])/deltaX[].n << endl;
	cout << "AbsY: " << (deltaY[]'*deltaY[])/deltaY[].n << endl;
	cout << "AbsZ: " << (deltaZ[]'*deltaZ[])/deltaZ[].n << endl;
	cout << "Rel: " << (Fx(r2,0,0)-Fx2(r2,0,0))/Fx2(r2,0,0) << endl;
	
	//v1x = Fx/rho*dt;
	//v1y = Fy/rho*dt;
	//v1z = Fz/rho*dt;
	
	v1x = v1x + Fx2/rho*dt;
	v1y = v1y + Fy2/rho*dt;
	v1z = v1z + Fz2/rho*dt;
	
	Ux = Ux +  v1x * dt;
	Uy = Uy +  v1y * dt;
	Uz = Uz +  v1z * dt;
	

	//just some animation stuff here
	//...
}
2

You have to be careful when switching from problem to varf, you have to compute the opposite of your linear form. Have you checked that? piezoTrans and piezoTrans2 seems to be exactly the same, so any code using varf (you can check without PETSc) will not work.

3

Sorry, I don’t really understand what you mean by “the opposite of my linear form”.

4

Multiply by -1. See FreeFem-sources/Laplace.edp at develop · FreeFem/FreeFem-sources · GitHub vs. FreeFem-sources/first.edp at develop · FreeFem/FreeFem-sources · GitHub.

5

Thank you very much! That finally solved my problem!