here is the code and the software told me there was an error.Please help me out!!
border lower(t=0,410){x=t;y=0;}
border right(t=0,91){x=410;y=t;}
border upper(t=410,0){x=t;y=91;}
border left(t=91,0){x=0;y=t;}
int c=6.7;//长度量纲单位均为mm
int solid=1;
int[int] YC=[1,3.3,5.5,7.7,9.9,12.1,14.3,16.5,18.7,20.9,23.1,25.3,27.5,29.7,31.9,34.1,36.3,38.5,40.7,42.9,45.1,47.3,49.5,51.7,53.9,56.1,58.3,60.5,62.7,64.9,67.1];
border s1(t=0,2pi;i)
{
x=YC[i]c+ccos(t);
y=YC[0]c+csin(t);
label=solid;
}
border s2(t=0,2pi;i)
{
x=YC[i]c+ccos(t);
y=YC[1]c+csin(t);
label=solid;
}
border s3(t=0,2pi;i)
{
x=YC[i]c+ccos(t);
y=YC[2]c+csin(t);
label=solid;
}
border s4(t=0,2pi;i)
{
x=YC[i]c+ccos(t);
y=YC[3]c+csin(t);
label=solid;
}
border s5(t=0,2pi;i)
{
x=YC[i]c+ccos(t);
y=YC[4]c+csin(t);
label=solid;
}
border s6(t=0,2pi;i)
{
x=YC[i]c+ccos(t);
y=YC[5]c+csin(t);
label=solid;
}
border s7(t=0,2*pi;i)
{
x=YC[i]c+ccos(t);
y=YC[6]c+csin(t);
label=solid;
}
int[int] nn=[100,20,100,20];
int[int] mm(31);
mm=-3;
int[int] ll(31);
ll=10;
mesh th1=buildmesh(lower(10)+right(10)+upper(10)+left(10)+s1(mm)+s2(mm)+s3(mm)+s4(mm)+s5(mm)+s6(mm)+s7(mm));
mesh th2= buildmesh(s1(ll)+s2(ll)+s3(ll)+s4(ll)+s5(ll)+s6(ll)+s7(ll));
mesh th=th1+th2;
plot(th,wait=1,ps=“total.eps”);
plot(th1,wait=1,ps=“pore.eps”);
plot (th2, wait=1,ps=“solid.eps”);
real s0=clock();
fespace v0h(th1,P1);
fespace v1h(th2,P1); //温度等变量空间
fespace v0h1(th1,P1b);
fespace v1h1(th2,P1b); //速度空间
fespace v0h2(th1,P0);
fespace v1h2(th2,P0);//梯度空间
v0h u1,v1,up1=0,up2=0;
v1h u2,v2,up3=0,up4=0;
macro def(u) [1/2*[2dx(u),dy(u)],1/2[dy(u),0]] //EOM
macro Grad(u) [dx(u),dy(u)] // EOM
macro UGrad(u1,u2) [u1*dx(u2)+u1*dy(u2)] //EOM
macro Mo(u) [pow(dx(u)*dx(u)+dy(u)*dy(u),1/2)] //EOM
v0h duf,uf1,uf=0.00001,vf1,vf=0.00001;
v1h dus,us=0,us1;
v0h dpre,pre1,pre=0.00001;
v0h dwf,wf1,wf=0.00001;
v0h dfi,fi1,fi=0.00001;
v0h setaf1=273.15;
v1h setas1=273.15;
//参数改单位
real alpha=0;
real rawl=1000; //kg/m3
real rawi=913;
real raws=2647;
real c1=1;
real dt=0;
func pfi=(6pi^2-15fi+10)fi^3;
func dpfi=30fi^4-60fi^3+30fi^2;
func gfi= pow(fi,2)pow(1-fi,2);
func dgfi=2fi+4fi^3-6fi^2;
func pfi1=(6fi1^2-15fi1+10)fi1^3;
func dpfi1=30fi1^4-60fi1^3+30fi1^2;
func gfi1= pow(fi1,2)pow(1-fi1,2);
func dgfi1=2fi1+4fi1^3-6fi1^2;
func rawf=pfirawl+(1-pfi)rawi;
real kl=2e10;
real ki=5.56e9;
real ks=5e10;// bulk modulus
real ass=10e-6;
real asi=5.1e-5;
real asl=1.8e-4; //heat dilation coefficient
real asf=pfiasl+(1-pfi)asi;
real mui=4.17e9;
real mus=37.5e9; //shear modulus
real muf=1;//viscosity
real seff=273.15;//k
real sefs=273.15;
real avi=1.09;// volume expansion due to ice formation
real D=1e-8;//水动力弥散系数
real nos=0.58;
real kos=1e-13;
real g=9.8;
matrix I=
[[1,0],[0,1]];
real hl=0.611; //conductivity
real hi=2.222;
real hs=7.694;
func hf=pfihl+(1-pfi)hi;
real cl=4200; //heat capacity
real ci=1900;
real cs=740;
func cf=pficl+(1-pfi)ci;
real sl=3.886;//标准熵
real si=2.6605;
func ns=nos(1-dx(us)-dy(us)+3ass*(u2-setas1));
func nf=1-ns;//后期可能需要变为循环语句
func M=1/((sl-si)rawf);
real rs=0.5;
func rf=pow(nf/ns,1/3)rs;
func W=3rf;
real ita=pow(6rf,1/2);
real kfss=0.55;
func kfsv=3nfkfss/rf;
func raws1=raws*(1-3ass(u2-setas1));
real a=0.003;
problem pref(pre,dpre)=int2d(th1)((dx(dpre)dx(pre)+dy(dpre)dy(pre))/(rawfc1^2)) +int2d(th1)(-vf1dx(dpre)-vf1dy(dpre))
+int1d(th1)(dpre(N.xvf1+N.yvf1))+on(lower,right,upper,left,pre=0);
problem momf(wf,dwf)=int2d(th1)(dx(dwf)((pfirawf/rawl)(kl(dx(uf1)+dy(uf1))-(3klasl*(setaf1-seff)))+((1-pfi)rawf/rawi)((2mui(dx(uf1)+dy(uf1))+ki*(dx(uf1)+dy(uf1))-3kiasi*(setaf1-seff))+nfpre))+
dy(dwf)((pfirawf/rawl)(kl*(dx(uf1)+dy(uf1))-(3klasl*(setaf1-seff)))+((1-pfi)rawf/rawi)((2mui(dx(uf1)+dy(uf1))+ki*(dx(uf1)+dy(uf1))-3kiasi*(setaf1-seff))+nfpre)))
-int2d(th1)(dwfrawfg)
-int2d(th1)(nfmufwfdwf/kos)
-int1d(th1)(dwf*(N.x*((pfirawf/rawl)(kl*(dx(uf1)+dy(uf1))-(3klasl*(setaf1-seff)))+((1-pfi)rawf/rawi)((2mui(dx(uf1)+dy(uf1))+ki*(dx(uf1)+dy(uf1))-3kiasi*(setaf1-seff))+nfpre))+
N.y((pfirawf/rawl)(kl*(dx(uf1)+dy(uf1))-(3klasl*(setaf1-seff)))+((1-pfi)rawf/rawi)((2mui(dx(uf1)+dy(uf1))+ki*(dx(uf1)+dy(uf1))-3kiasi*(setaf1-seff))+nfpre))))
+on(lower,right,upper,left,wf=0); //蓝色部分为液相有效应力
problem velof(uf,duf)=int2d(th1)(2dufalphauf)+int2d(th1)(duf*(wf1dx(uf)+wf1dy(uf)))-int2d(th1)(dufalphaconvect([up1,up2],-dt,wf1))+on(lower,right,upper,left,uf=0.0001);
problem moms(us,dus)=int2d(th)(dus*((pfirawf/rawl)(kl*(dx(uf1)+dy(uf1))-(3klasl*(setaf1-seff)))+((1-pfi)rawf/rawi)(2mui(dx(uf1)+dy(uf1))+ki*(dx(uf1)+dy(uf1))-3kiasi*(setaf1-seff))))
+int2d(th)(dusraws(2mus(dx(us)+dy(us))+ks*(dx(us)+dy(us))))
-int2d(th)(dus*(3ksass*(setas1-sefs)-(nssetas1/setaf1-ns+1)pre1))
+int2d(th)(dus(raws+rawf)g)
-int1d(th)(dus(N.x((pfirawf/rawl)(kl*(dx(uf1)+dy(uf1))-(3klasl*(setaf1-seff)))+((1-pfi)rawf/rawi)(2mui(dx(uf1)+dy(uf1))+ki*(dx(uf1)+dy(uf1))-3kiasi*(setaf1-seff)))))
-int1d(th)(dusN.x(2mus(dx(us)+dy(us))+ks*(dx(us)+dy(us))))
+int1d(th)(dusN.x(3ksass*(setas1-sefs)-nspre1setas1/setaf1+nspre1-pre1))
-int1d(th)(dus(N.y*((pfirawf/rawl)(kl*(dx(uf1)+dy(uf1))-(3klasl*(setaf1-seff)))+((1-pfi)rawf/rawi)(2mui(dx(uf1)+dy(uf1))+ki*(dx(uf1)+dy(uf1))-3kiasi*(setaf1-seff)))))
-int1d(th)(dusN.y(2mus(dx(us)+dy(us))+ks*(dx(us)+dy(us))))
+int1d(th)(dusN.y(3ksass*(setas1-sefs)-nspre1setas1/setaf1+ns*pre1-pre1))
+on(lower,right,upper,left,us=0.0000001);
problem order(fi,dfi)=int2d(th1)(dfialphafi/M)
-int2d(th1)(alphadficonvect([up1,up2],-dt,fi1)/M)
//-int2d(th1)(dfidpfi(rawf*(cl-ci)(setaf1log(setaf1/273.15)-setaf1+seff))
//+int2d(th1)(dfidpfirawf*(sl-si)seff(setaf1-seff)/seff)
//+int2d(th1)(0.5dfidpfiitaita*(pow(dy(fi),2)+pow(dx(fi),2))(rawi-rawl)/rawf)
//-int2d(th1)(dfigfidpfiW*(rawi-rawl)/rawf)
//+int2d(th1)(dfidgfidpfiW)
//+int2d(th1)(dfidpfipre1(rawl-rawi))
+int2d(th1)(itaita(dx(fi)dx(dfi)+dy(fi)dy(dfi)))
+int1d(th1)(N.xdfiitaita(dx(fi)+dy(fi)))
+int1d(th1)(N.ydfiitaita(dx(fi)+dy(fi)))
+on(lower,right,upper,left,fi=1);// fi=1,与非线性相关的项为与dpfi相关的所有项
problem temps(u2,v2)=
int2d(th2)(alphav2(nspreu2/setaf1)(dx(us)+dy(us)))
-int2d(th2)(alphav2*(nspreu2/setaf1)convect([up3,up4],-dt,(dx(us1)+dy(us1))))
+int2d(th2)(alphav23ksass(dx(us)+dy(us)))
-int2d(th2)(alphav23ksassconvect([up3,up4],-dt,(dx(us1)+dy(us1))))
+int2d(th2)(alphav2rawscsu2)
-int2d(th2)(alphav2rawscsconvect([up3,up4],-dt,setas1))
//+int2d(th2)(alphav26ns*(raws1/raws)(u2/setaf1)preassu2)
//+int2d(th2)(alphav26ns(raws1/raws)(u2/setaf1)preassconvect([up3,up4],-dt,setas1))
//-int2d(th2)(alphav218nspow(u2,2)prepow(ass,2)u2/setaf1)
//+int2d(th2)(alphav218nspow(u2,2)prepow(ass,2)convect([up3,up4],-dt,setaf1)/setaf1)
//-int2d(th2)(alphav2pow(u2,2)pow(setaf1,-2)pre3assnssetaf1)
//+int2d(th2)(alphav2pow(u2,2)pow(setaf1,-2)pre3assnsconvect([up3,up4],-dt,setaf1))
//-int2d(th2)(alphav2pow(u2,2)pow(setaf1,-1)3assnspre)
//+int2d(th2)(alphav2pow(u2,2)pow(setaf1,-1)3assnsconvect([up3,up4],-dt,pre1))
//+int2d(th2)(v2kfsv(u2-setaf1))
+int2d(th2)(v2nfmufwfwf/kos)
+int2d(th2)(hs*(dx(u2)dx(v2)+dy(u2)dy(v2)))
-int1d(th2)(N.xv2(hs*(dx(u2)+dy(u2))))
-int1d(th2)(N.yv2(hs*(dx(u2)+dy(u2))))
+int2d(th2,left,right)(v2N.xu2)
+int2d(th2,upper)(v2N.yu2);
+int2d(th2,lower)(v2*(a*(263.15-u2)-100));
solve tempf(u1,v1)=int2d(th1)(alphav1rawfcfu1)+int2d(th1)(alphav1rawfcfconvect([up1,up2],-dt,setaf1))
//+int2d(th1)(alphav1nfpre9asfasfu1u1)+int2d(th1)(alphav1nfpre9asfasfu1convect([up1,up2],-dt,setaf1))
+int2d(th1)(alphav1preasfu1)+int2d(th1)(alphav1preasfconvect([up1,up2],-dt,setaf1))
-int2d(th1)(alphav13nfasfu1pre)+int2d(th1)(alphav13nfasfu1convect([up1,up2],-dt,pre1))
+int2d(th1)(alphav1nfpre(dx(uf)+dy(uf)))
-int2d(th1)(alphav1nfpreconvect([up1,up2],-dt,(dx(uf1)+dy(uf1))))
+int2d(th1)(alphav1rawf*(pfiu13klasl/rawl+(1-pfi)u13kiasi/rawi)(dx(uf)+dy(uf)))
-int2d(th1)(alphav1rawf(pfiu13klasl/rawl+(1-pfi)u13kiasi/rawi)convect([up1,up2],-dt,(dx(uf1)+dy(uf1))))
+int2d(th1)(alphav1Wdgfifi)
-int2d(th1)(alphav1Wdgficonvect([up1,up2],-dt,fi1))
+int2d(th1)(v1itaita(dx(fi)dx(dfi)+dy(fi)dy(dfi)))
+int2d(th1)(alphav1(rawl-rawi)0.5itaita(dx(fi)dx(fi)+dy(fi)dy(fi))dpfinfrawffi/pow(rawf,2))
-int2d(th1)(alphav1(rawl-rawi)0.5itaita(dx(fi)dx(fi)+dy(fi)dy(fi))dpfinfrawfconvect([up1,up2],-dt,fi1)/pow(rawf,2))
+int2d(th1)(alphav1(rawl-rawi)Wgfidpfinfrawffi/pow(rawf,2))
-int2d(th1)(alphav1(rawl-rawi)Wgfidpfinfrawfconvect([up1,up2],-dt,fi1)/pow(rawf,2))
+int2d(th1)(alphav1(rawl-rawi)Wgfidpfinfrawffi/pow(rawf,2))
-int2d(th1)(alphav1(rawl-rawi)Wgfidpfinfrawfconvect([up1,up2],-dt,fi1)/pow(rawf,2))
+int2d(th1)(-alphav13preu1*(aslrawi-asirawl)dpfinfrawffi/pow(rawf,2))
-int2d(th1)(-alphav13preu1*(aslrawi-asirawl)dpfinfrawfconvect([up1,up2],-dt,fi1)/pow(rawf,2));
+int2d(th1)(-alphav1muipow(dx(uf)+dy(uf),2)dpfinfrawffi/rawi)
-int2d(th1)(-alphav1muipow(dx(uf)+dy(uf),2)dpfinfrawfconvect([up1,up2],-dt,fi1)/rawi)
+int2d(th1)(alphav1pow(dx(uf)+dy(uf),2)0.5(kl/rawl-ki/rawi)dpfinfrawffi)
-int2d(th1)(alphav1pow(dx(uf)+dy(uf),2)0.5(kl/rawl-ki/rawi)dpfinfrawfconvect([up1,up2],-dt,fi1)/rawi)
+int2d(th1)(alphav1(dx(uf)+dy(uf))seff3*(aslkl/rawl-asiki/rawi)dpfinfrawffi)
-int2d(th1)(alphav1(dx(uf)+dy(uf))seff3*(aslkl/rawl-asiki/rawi)dpfinfrawfconvect([up1,up2],-dt,fi1)/rawi)
+int2d(th1)(alphav1(dx(uf)+dy(uf))seff3*(aslkl/rawl-asiki/rawi)dpfinfrawffi)
-int2d(th1)(alphav1(dx(uf)+dy(uf))seff3*(aslkl/rawl-asiki/rawi)dpfinfrawfconvect([up1,up2],-dt,fi1)/rawi)
+int2d(th1)(alphav1(cl-ci)(u1-seff)dpfinfrawffi)
-int2d(th1)(alphav1*(cl-ci)(u1-seff)dpfinfrawfconvect([up1,up2],-dt,fi1)/rawi)
+int2d(th1)(alphav1*(sl-si)seffdpfinfrawffi)
-int2d(th1)(alphav1*(sl-si)seffdpfinfrawfconvect([up1,up2],-dt,fi1)/rawi)
+int2d(th1)(alphav13kiavi(dx(uf)+dy(uf))pow(rawi,-1)dpfinfrawffi)
-int2d(th1)(alphav13kiavi(dx(uf)+dy(uf))pow(rawi,-1)dpfinfrawfconvect([up1,up2],-dt,fi1)/rawi)
-int2d(th1)(v1kfsv*(u1-u2))
+int2d(th1)(hs*(dx(u1)dx(v1)+dy(u1)dy(v1)))
-int1d(th1)(N.xv1(hs*(dx(u1)+dy(u1))))
-int1d(th1)(N.yv1(hs*(dx(u1)+dy(u1))))
+int2d(th1,left,right)(v1N.xu1)
+int2d(th1,upper)(v1N.yu1)
+int2d(th1,lower)(v1*(a*(263.15-u1)-100));
//+on(left,right,N.xu1=0)
//+on(lower,N.yu1=0)
//+on(upper,0.003(263.15-u1)=100);
plot (coef=0.02,cmm="u2 et pre ",pre,u2,ps=“test.eps”,value=1,wait=1);
{
real[int] tt(21),pp(21),yy(21);
for (int i=0;i<21;i++)
{
yy[i]=i/20;
tt[i]=u2(200,i/20.);
pp[i]=pre(i/20.,45);
}
cout << “” << yy << endl;
}
dt=0.5;
for (int i=0;i<=500;i++)
{
alpha=1/dt;
for(int j=0;j<=500;j++)
{
up1=vf;
up2=0;
up3=vf-wf;
up4=0;
pre1=pre;
wf1=wf;
uf1=uf;
us1=us;
fi1=fi;
setaf1=u1;
setas1=u2;
pref;plot(pre,wait=1);
momf;plot(wf,wait=1);
velof;plot(vf,wait=1);
moms;plot(us,wait=1);
order;plot(fi,wait=1);
temps;
tempf;
}
}