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
//...
}
```