///////////////////////////////// // IPOPT3 load "ff-Ipopt"; mesh Th=square(50,50); fespace Xh(Th,[P2,P2,P1]); fespace Vh(Th,P2), Qh(Th,P1); Vh lbu=-1.e19; Vh lbv=-1.e19; Qh lbp=-1.e19; Vh ubu=1.e19; Vh ubv=1.e19; Qh ubp=1.e19; Vh u,v,uu,vv,du,dv,uk,vk, uold=0,vold=0; Qh p,pp,dp,pk; real nu=0.1, T=1., dt = 0.1; int m, M= T/dt; func real J(real[int]& Y) { Xh [u,v,p]; u[]=Y; // original J, DJ is the correct variaitional of J // using this J form, the line search function is ON for the IPOPT below return int2d(Th)( 0.5*(u^2 + v^2)/dt + 0.5*nu*(dx(u)*dx(u) + dy(u)*dy(u) + dx(v)*dx(v) + dy(v)*dy(v)) - 0.5*(p^2)*1.e-1 - (uold*u+vold*v)/dt); // // modified J, DJ is NOT the correct variaitional of J // // using thia J form, the line search function MUST turn OFF for the IPOPT below // return int2d(Th)( 0.5*(u^2 + v^2)/dt + 0.5*nu*(dx(u)*dx(u) + dy(u)*dy(u) // )); } func real[int] DJ(real[int]& Y) { Xh [u,v,p]; u[]=Y; varf aa ([uk,vk,pk],[uu,vv,pp]) = int2d(Th)( (u*uu+v*vv)/dt + nu*(dx(u)*dx(uu) + dy(u)*dy(uu) + dx(v)*dx(vv) + dy(v)*dy(vv)) - p*pp*1.e-1 ) - int2d(Th)((uold*uu+vold*vv)/dt); real[int] G= aa(0,Xh); return G;} matrix H;//global vairable for Hessien matrix Overwise => seg fault in Ipopt func matrix HJ(real[int]& Y) { Xh [u,v,p]; u[]=Y; varf AK ([uu,vv,pp],[du,dv,dp]) = int2d(Th)( (du*uu+dv*vv)/dt + nu*(dx(du)*dx(uu) + dy(du)*dy(uu) + dx(dv)*dx(vv) + dy(dv)*dy(vv)) - dp*pp*1.e-1 ); H = AK(Xh,Xh); return H;} varf vGamma1(u,uu) = on(1,2,4,u=1); real[int] onGamma1 = vGamma1(0,Vh); lbu[] = onGamma1 ? 0.:lbu[]; // for u lbv[] = onGamma1 ? 0.:lbv[]; // for v ubu[] = onGamma1 ? 0.:ubu[]; // for u ubv[] = onGamma1 ? 0.:ubv[]; // for v varf vGamma3(u,uu) = on(3,u=1); real[int] onGamma3 = vGamma3(0,Vh); lbu[] = onGamma3 ? 1.:lbu[]; // for u lbv[] = onGamma3 ? 2.7:lbv[]; // for v ubu[] = onGamma3 ? 1.:ubu[]; // for u ubv[] = onGamma3 ? 2.7:ubv[]; // for v Xh [lb1,lb2,lb3]=[lbu,lbv,lbp]; Xh [ub1,ub2,ub3]=[ubu,ubv,ubp]; Xh [w1,w2,wp]=[0,0,0]; for(m=0;m