Dear @prj
I have taken a simple model for Navier Stoke Equation over a square.
which I have taken from Navier-Stoke.
- I am getting an error as " operation error ".
code is:
real nu = 1.0/50.;
real eps = 1e-6;
verbosity = 0;
mesh Th = square(10,10);
fespace Xh(Th, P2);
Xh u1, u2;
Xh du1, du2;
Xh v1, v2;
fespace Mh(Th, P1);
Mh p;
Mh q;
Mh dp;
macro Grad(u1, u2) [ dx(u1), dy(u1), dx(u2), dy(u2)] //
macro UgradV(u1, u2, v1, v2) [ [u1, u2 ]' * [ dx(v1), dy(v1) ], [u1, u2]' * [dx(v2), dy( v2)] ] //
macro div( u1, u2) ( dx(u1) + dy(u2) ) //
// initialization
u1 = x^2 + y^2 ;
u2 = 0;
// viscosity loop
while(1)
{
int n;
real err = 0 ;
for ( n =0; n< 10; n++)
{
solve viscous([ du1, du2, dp ], [ v1, v2, q] ) =
int2d(Th) ( nu *( Grad(u1, u2)' * Grad( v1, v2))
- div( v1, v2)*dp
- div( du1, du2)*q
+ UgradV( du1, du2, u1, u2 )'* [v1, v2]
+ UgradV( u1, u2, du1, du2 )'* [ v1, v2 ]
- 1e-8*dp*q )
- int2d( Th) ( nu * (Grad(u1, u2)'Grad(v1, v2))
+ UgradV( u1, u2, u1, u2)' * [ v1, v2 ]
- p * div(v1, v2 ) - q * div(u1, u2)
)
+ on( 1,2,4, du1 = 0 , du2 = 0 )
+ on ( 3, du1 = 1, du2 = 0 )
u1[] -= du1[];
u2[] -= du2[];
p[] -= dp[];
real Lu1=u1[].linfty, Lu2=u2[].linfty, Lp=p[].linfty;
err = du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
cout << n << " err = " << err << " " << eps << " rey = " << 1./nu << endl;
if(err < eps) break; //converge
if( n>3 && err > 10.) break; //blowup
}
}
Kindly, suggest me to fix it.