Hello, teachers. I’m running the thermally coupled MHD equations and I’m having the following problem. The display is out of memory, when I reduce the mesh size it works but the mesh size is not very large. Please tell me why. Here is codes
load "UMFPACK64"
//load "PETSc_FreeFem"
//定义微元
load "splitmesh3"
load "Element_P3"
macro div(u1,u2) ( dx(u1) + dy(u2) ) //
macro dot(u1,u2,v1,v2) ( (u1) * (v1) + (u2) * (v2) ) //
macro Ugradv1(u1,u2,v1) ( (u1) * dx(v1) + (u2) * dy(v1)) //
macro trilinear(u1,u2,v1,v2,w1,w2) ( dot(Ugradv1(u1,u2,v1),Ugradv1(u1,u2,v2),w1,w2)) //
macro Bu(u1,u2,v1,v2) ( dx(u1) * dx(v1) + dx(u2) * dx(v2) +dy(u1) * dy(v1) + dy(u2) * dy(v2) ) //B_u
macro mo1(u1,u2,v1,v2,w1,w2) (dot(u1,u2,v1,v2)*div(w1,w2))//
func real f1(real t) {return 0;}
func real f2(real t) {return 0;}
func real g1(real t) {return 0;}
func real g2(real t) {return 0;}
func real Psi(real t) {return 0;}
int i=0;
//给出参数
real nI,nu=0,mu=0,alpha1=0.3,alpha2=0,alpha3=0,s=1.0,kappa=1.0,beta1=0,beta2=-1.;//有Voigt
real t;
real [int] arr=[64];
//空间循环
for(i=0;i<1;i++){
nI=arr[i];
//建立网格
mesh ThF=square(nI,nI,[-pi+2.*pi*x,-pi+2.*pi*y]);
ThF = splitmesh3(ThF);
int freedom = ThF.nv;
cout<<"\n"<<"freedom="<<freedom;
plot(ThF,wait=true);
//建立空间
fespace Xh(ThF,[P2,P2,P1dc,P2,P2,P1dc,P2],periodic=[[4,y],[2,y],[1,x],[3,x]]);
fespace MhI(ThF,P2,periodic=[[4,y],[2,y],[1,x],[3,x]]);//温度
MhI curlB;
Xh [u1h,u2h,phmh,b1h,b2h,qhmh,thetah],[u1ho,u2ho,phmho,b1ho,b2ho,qhmho,thetaho],[u1hoo,u2hoo,phmhoo,b1hoo,b2hoo,qhmhoo,thetahoo],[v1h,v2h,lambdah,chi1h,chi2h,rh,phih];
//MhI u1ho,u2ho,u1hoo,u2hoo,b1ho,b2ho,b1hoo,b2hoo,thetaho,,dxb2ho;
int freedom1 = Xh.ndof;
int freedom2 = MhI.ndof;
cout<<"\n"<<"freedom1="<<freedom1;
cout<<"\n"<<"freedom2="<<freedom2;
//给定初值
real t;
t=0.0;
[u1ho,u2ho,phmho,b1ho,b2ho,qhmho,thetaho]=[-sin(y+2.),sin(x+1.4),0, -(1./3.)*sin(y+6.2), (2./3.)*sin(2.*x+2.3),0,0];
[u1hoo,u2hoo,phmhoo,b1hoo,b2hoo,qhmhoo,thetahoo]=[-sin(y+2.),sin(x+1.4),0, -(1./3.)*sin(y+6.2), (2./3.)*sin(2.*x+2.3),0,0];
//thetaho = 0;
real tmmh,tmp1;
int m;
real TIME=2.7;
real deltat=0.01;
real maxstepnum=TIME/deltat;
//时间循环
for(m=1;m<= maxstepnum;m++){
cout<<"\n"<<"#######################";
cout<<"\n"<<"m="<<m;
t=deltat*m;
tmmh=deltat*(m-0.5);
cout<<"\n"<<"t="<<t;
solve TMHD([u1h,u2h,phmh,b1h,b2h,qhmh,thetah],[v1h,v2h,lambdah,chi1h,chi2h,rh,phih],solver=UMFPACK)=
int2d(ThF)(
(1.0/deltat)*dot(u1h,u2h,v1h,v2h)
+(0.5*nu+alpha1^2/deltat)*Bu(u1h,u2h,v1h,v2h)
+0.5*trilinear(1.5*u1ho-0.5*u1hoo,1.5*u2ho-0.5*u2hoo,u1h,u2h,v1h,v2h)
-s*0.5*trilinear(1.5*b1ho-0.5*b1hoo,1.5*b2ho-0.5*b2hoo,b1h,b2h,v1h,v2h)
-phmh*div(v1h,v2h)
-s*0.25*dot(1.5*b1ho-0.5*b1hoo,1.5*b2ho-0.5*b2hoo,b1h,b2h)*div(v1h,v2h)
-0.5*thetah*dot(beta1,beta2,v1h,v2h)
+0.000001*phmh*lambdah
)
+int2d(ThF)(
-(1.0/deltat)*dot(u1ho,u2ho,v1h,v2h)
+(0.5*nu-alpha1^2/deltat)*Bu(u1ho,u2ho,v1h,v2h)
+0.5*trilinear(1.5*u1ho-0.5*u1hoo,1.5*u2ho-0.5*u2hoo,u1ho,u2ho,v1h,v2h)
-s*0.5*trilinear(1.5*b1ho-0.5*b1hoo,1.5*b2ho-0.5*b2hoo,b1ho,b2ho,v1h,v2h)
-s*0.25*dot(b1ho,b2ho,b1ho,b2ho)*div(v1h,v2h)
-0.5*thetaho*dot(beta1,beta2,v1h,v2h)
-dot(f1(tmmh),f2(tmmh),v1h,v2h)
)
+int2d(ThF)(
div(u1h,u2h)*lambdah
)
+int2d(ThF)(
(1.0/deltat)*dot(b1h,b2h,chi1h,chi2h)
+(0.5*mu+alpha2^2/deltat)*Bu(b1h,b2h,chi1h,chi2h)
-0.5*trilinear(1.5*b1ho-0.5*b1hoo,1.5*b2ho-0.5*b2hoo,u1h,u2h,chi1h,chi2h)
+0.5*trilinear(1.5*u1ho-0.5*u1hoo,1.5*u2ho-0.5*u2hoo,b1h,b2h,chi1h,chi2h)
-qhmh*div(chi1h,chi2h)
+0.000001*qhmh*rh
)
+int2d(ThF)(
-(1.0/deltat)*dot(b1ho,b2ho,chi1h,chi2h)
+(0.5*mu-alpha2^2/deltat)*Bu(b1ho,b2ho,chi1h,chi2h)
-0.5*trilinear(1.5*b1ho-0.5*b1hoo,1.5*b2ho-0.5*b2hoo,u1ho,u2ho,chi1h,chi2h)
+0.5*trilinear(1.5*u1ho-0.5*u1hoo,1.5*u2ho-0.5*u2hoo,b1ho,b2ho,chi1h,chi2h)
-dot(g1(tmmh),g2(tmmh),chi1h,chi2h)
)
+int2d(ThF)(
div(b1h,b2h)*rh
)
+int2d(ThF)(
(1.0/deltat)*thetah*phih
+(0.5*kappa+alpha3^2/deltat)*dot(dx(thetah),dy(thetah),dx(phih),dy(phih))
+0.5*dot(1.5*u1ho-0.5*u1hoo,1.5*u2ho-0.5*u2hoo,dx(thetah),dy(thetah))*phih
)
+int2d(ThF)(
-(1.0/deltat)*thetaho*phih
+(0.5*kappa-alpha3^2/deltat)*dot(dx(thetaho),dy(thetaho),dx(phih),dy(phih))
+0.5*dot(1.5*u1ho-0.5*u1hoo,1.5*u2ho-0.5*u2hoo,dx(thetaho),dy(thetaho))*phih
-Psi(tmmh)*phih
);
[u1hoo,u2hoo,phmhoo,b1hoo,b2hoo,qhmhoo,thetahoo]=[u1ho,u2ho,phmho,b1ho,b2ho,qhmho,thetaho];
[u1ho,u2ho,phmho,b1ho,b2ho,qhmho,thetaho]=[u1h,u2h,phmh,b1h,b2h,qhmh,thetah];
curlB=dx(b2h)-dy(b1h);
plot(curlB);
cout << " *** real time = " << clock() << endl;
}
cout<<"t="<<t<<"\n";
int i,j;
ofstream ffB("Orszag-Tang_curlB_48_0.3_sv.dat");
ffB<<"Variables="<<"X"<<","<<"Y"<<","<<"curlB"<<endl;
ffB<<"Zone"<<" "<<"N="<<ThF.nv<<","<<"E="<<ThF.nt<<","<<"F=FEPOINT,ET=TRIANGLE"<<endl;
for(i=0;i<ThF.nv;i++)
{
ffB<<ThF(i).x<<" "<<ThF(i).y<<" "<<curlB(ThF(i).x,ThF(i).y)<<endl;
}
for(i=0;i<ThF.nt;i++)
{
for(j=0;j<3;j++)
{
ffB<<ThF[i][j]+1<<" ";
}
if(j==3)
{
ffB<<endl;
}
}
}
Here is output
-- Square mesh : nb vertices =4225 , nb triangles = 8192 , nb boundary edges 256
freedom=12417
freedom1=393216
freedom2=49152
#######################
m=1
t=0.01 Error Umfpack -1 : out_of_memory current line = 122
Exec error : Error Umfpack -1 : out_of_memory
-- number :1
catch an erreur in solve => set sol = 0 !!!!!!!
Exec error : Error Umfpack -1 : out_of_memory
-- number :1
err code 8 , mpirank 0
Thank you!!!