Define a convection constant vectro!

General Discussion
1

Hi everybody,
If you allow, I have a question about defining a constant advection vector for my reaction-advection-diffusion system. More precisely, I have the following program where the vector director of the direction of advection is given by the gradient of “c”. My question is to define
a constant vector [v1,v2] in sorts that the direction is from outside the circle to the part of boundary with the angle between (0,theta). Thank you very much in advance. The program is next:
real Lx = 3;
real Ly = 2;
real dt = 0.05;
real uf = 1;
real rhoc = 100;
// Define mesh boundary
real theta = pi/10.;
//real theta2 = 5.pi/11;
border C(t=0. , theta){x=cos(t); y=sin(t);}
border D(t=theta , 2.pi){x=cos(t); y=sin(t);}
// The triangulated domain Th is on the left side of its boundary
mesh Th = buildmesh(C(150)+D(100));
//mesh Th = square(60, 40);
plot(Th);//, ps = “mesh.eps”);
fespace Uh(Th, P2);
fespace Vh(Th, P1);
Uh rho, lrho, rhoold, rhoh, c, ch, cold;
Vh u1, u2, u, n1, n2, v1, v2, rhop1, cp1;
rho = rhocexp( -40(x-Lx/4)^2 - 40*(y-Ly/2)^2 )

  • rhocexp( -100(x-3Lx/4)^2 - 40(y-Ly/2)^2 )
  • rhocexp( -100(x-0.55Lx)^2 - 40(y-Ly/2)^2 );
    Th = adaptmesh(Th, rho);
    //Th = adaptmesh(Th, rho, periodic=[[3,x],[1,x]]);
    plot(Th);
    rho = rho;
    c = rho/rhoc;
    rhoold = rho;
    cold = c;
    // plot (rho, nbiso=50, wait=1);
    //
    problem migr([rho, c], [rhoh, ch]) =
    int2d(Th)(rhorhoh/dt)
    -int2d(Th)(convect([v1,v2], -dt,rhoold)rhoh/dt)
    +int2d(Th)(dx(v1)rhorhoh+dy(v2)rhorhoh)
    +int2d(Th)(0.01    dx(rho)dx(rhoh)+0.01dy(rho)dy(rhoh))
    -int2d(Th)(0.01    rho    (rhoc-rhoold)rhoh)
    +int2d(Th)(c    ch/dt)
    -int2d(Th)(cold*ch/dt)
    +int2d(Th)(dx(c)*dx(ch)+dy(c)dy(ch))
    -int2d(Th)(10    (rho/rhoc-c)*ch);
    //
    for (int it=0; it<20; it++) {
    for (int substep=0; substep<2; substep++){
    //u1 = -dx(cold);
    //u2 = -dy(cold);
    u1 = dx(cold);
    u2 = dy(cold);
    v1 = 0.5 * u1 * rhoold/rhoc;
    v2 = 0.5 * u2 * rhoold/rhoc;
    migr;
    Th = adaptmesh(Th, rho);
    rho=rho;
    c=c;
    rhoold = rho;
    cold = c;
    }
    // Visu
    rhop1 = rho; cp1 = c;
    plot(Th, rhop1, nbiso=50, fill=0, value=1,
    ps=“migr_it=”+it+".eps");
    }
    cout << “Done.\n”;