Hello! I am not sure how to solve a block system with EPSSolve.
My left hand side is:
J= [ [AA, AB, 0 ],
[BA', BB, BC ],
[0 , BC', 0 ] ];
and my right hand side is
Mat<complex> Loc = [[M1M ,0 ,0 ],
[0 ,M2M ,0 ],
[0 ,0 ,M3M ]];
the related vector is T,[u1,u2],p for the blocks.
Could you please check out my script to solve a conjugate heat transfer eigenvalue problem:
// run with MPI: ff-mpirun -np 4 script.edp
// NBPROC 4
//load "iovtk"
load "PETSc-complex"
macro dimension()2//
include "macro_ddm.idp"
func PkU = [P2, P2];
func Pkp = P1;
func PkT = P2;
macro grad(u)[dx(u), dy(u)]//
macro div(u)(dx(u#1) + dy(u#2))//
macro UgradV(u, v)[[u#1, u#2]' * [dx(v#1), dy(v#1)],
[u#1, u#2]' * [dx(v#2), dy(v#2)]]//
int inner = 1;
int outer = 2;
border c1(t=0,2*pi){
x=0.5*cos(t);
y=0.5*sin(t);
label=inner;
}
border c2(t=0,2*pi){
x=1*cos(t);
y=1*sin(t);
label=outer;
}
mesh Th1=buildmesh(c1(-100)+c2(100));
mesh Th2=buildmesh(c1(100));
mesh Th3=Th1+Th2;
plot(Th3,wait=1);
fespace VhU(Th1,[P2,P2]);
fespace Vhp(Th1,P1);
fespace VhT(Th3,P2);
VhU<complex> [uc1,uc2]=[0,0];
VhT<complex> Tc=y+1;
plot(Th3,Tc,wait=1);
real Pr = 1;
real Ra = 1e3;
real tM=getARGV("-tM", 0.5); // thetaM
real kappa,nu;
func real DefCoef()
{
nu = sqrt(Pr/Ra);
kappa = 1./sqrt(Ra*Pr);
}
DefCoef();
varf A11(T,TT) = int2d(Th3)(
kappa * (grad(T)' * grad(TT))
+ [uc1,uc2]'*[dx(T), dy(T)]*TT)
+ on(outer, T=-y);
varf A12([u1,u2],[TT]) = int2d(Th1)(
[u1,u2]'*[dx(Tc), dy(Tc)]*TT)
//+ on(inner,outer, u1 = 0, u2 = 0)
;
varf A21([v1,v2],[T]) = int2d(Th1)(
- (T)*v2)
/* + on(outer, T=-y) */;
varf A22([u1,u2],[v1,v2]) = int2d(Th1)(
(UgradV(uc, u) + UgradV(u, uc))' * [v1, v2]
+ nu * (grad(u1)' * grad(v1) +
grad(u2)' * grad(v2)))
+ on(inner,outer, u1 = 0, u2 = 0);
varf A23([p],[v1,v2]) = int2d(Th1)(
- p * div(v))
;
varf mass1(T,TT) = int2d(Th3)( T*TT );
varf mass2([u1,u2],[v1,v2]) = int2d(Th1)( [u1,u2]'*[v1,v2] );
varf mass3(p,q) = int2d(Th1)( 0*p*q );
matrix<complex> A11m=A11(VhT,VhT);
matrix<complex> A12m=A12(VhU,VhT);
matrix<complex> A21m=A21(VhT,VhU);
matrix<complex> A22m=A22(VhU,VhU);
matrix<complex> A23m=A23(Vhp,VhU);
//Mat<complex> ;
Mat<complex> AA=A11m;
Mat<complex> AB=A12m;
Mat<complex> BA=A21m;
Mat<complex> BB=A22m;
Mat<complex> BC=A23m;
Mat<complex> J= [ [AA, AB, 0 ],
[BA', BB, BC ],
[0 , BC', 0 ] ];
matrix<complex> M1m = mass1(VhT, VhT);
matrix<complex> M2m = mass2(VhU, VhU);
matrix<complex> M3m = mass3(Vhp, Vhp);
Mat<complex> M1M = M1m;
Mat<complex> M2M = M2m;
Mat<complex> M3M = M3m;
Mat<complex> Loc = [[M1M ,0 ,0 ],
[0 ,M2M ,0 ],
[0 ,0 ,M3M ]];
Mat<complex> M(J, Loc, clean = true);
//fespace Wh(Th3, [P2, P2, P2, P1]); // complete space [u, v, T, p]
int dofGlobal = VhT.ndof + VhU.ndof + Vhp.ndof;
int nev = 10;
complex[int,int] vec(dofGlobal,nev);
complex[int] val(nev); // array to store eigenvalues
//Wh<complex>[int] [vec,vecB,vecC,vecD](nev); // array to store eigenvectors
complex s = getARGV("-shift_real", 1.0e-6) + getARGV("-shift_imag", 0.6) * 1i;
string params = "-eps_tol 1.0e-6 -eps_nev " + nev + " " +
"-eps_type krylovschur -eps_monitor_all " +
"-eps_target " + real(s) + "+" + imag(s) + "i";
int k = EPSSolve(J, M, array = vec, values = val, sparams = params);
I am not sure what’s wrong with my code.
For now, I will go back to use block system.