Convect function failing

General Discussion
1

Hi everyone, I am currently using FreeFEM to model the kuramoto sivashinsky equation
u_{t}+\nu u_{xxxx}+u_{xx}+uu_{x}=0
with periuodic spatial and temporal boundary conditions. I am modeling this equation by differentiating u once and then modeling both the ux and uy derivatives. I am trying to use the convect function to handle the u_{t} and uu_{x} terms however I always arrive at the error message:

error operator ( <7E_Array>, , <10LinearCombI7MGauche4C_F0E>…

I am not sure what I’m doing wrong here and would appreciate any advice that could be given. The code I am working on is listed below. Thanks in advance!

\begin{verbatim}
real nu = 1;
real lambda = 25;
real dt = 0.05;
real epsu = 1e-6;
//mesh
mesh Th = square(50,50);

//fespace

fespace Vh(Th, P2b,
periodic = [[0,lambda],[0,lambda]]);
Vh u=sin(lambdapix)+sin(lambdapiy);
Vh ux = dx(u);
Vh uy = dy(u);
Vh v;

for(int n = 0; n<100; n++){
Vh uxold = ux;
Vh uyold = uy;

//problem
solve ksx(ux, v, init = n, solver = CG, eps = epsu)
= int2d(Th) (uxv)
//+ int2d(Th) (
//convect([dx(ux), dy(ux)], -dt, ux)v) //convect term causes crash
+ int2d(Th)(
dtnu/2(dxx(ux)*dxx(v)+dyy(ux)dyy(v)))
- int2d(Th)(
dt/2(dx(uxold)*dx(v)+dy(uxold)dy(v)))
- int2d(Th)(
dt/2(dx(ux)dx(v)+dy(ux)dy(v)))
+ int2d(Th)(
dtnu/2(dxx(uxold)*dxx(v)+dyy(uxold)*dyy(v)))

		;
solve ksy(uy, v, init = n, solver = CG, eps = epsu)
  		= int2d(Th) (uy*v)
  		//+ int2d(Th) (
          		// convect([dx(uy), dy(uy)], -dt, uy)*v) //convect term causes crash
  		+ int2d(Th)(
         		 dt*nu/2*(dxx(uy)*dxx(v)+dyy(uy)*dyy(v)))
  		- int2d(Th)(
          		dt/2*(dx(uyold)*dx(v)+dy(uyold)*dy(v)))
 		- int2d(Th)(
          		dt/2*(dx(uy)*dx(v)+dy(uy)*dy(v)))
  		+ int2d(Th)(
        		 dt*nu/2*(dxx(uyold)*dxx(v)+dyy(uyold)*dyy(v)))

                   ;

        //solver here should be laplace u = gradient f where f is (ux,uy), after this I can plot u
        solve ksu(u,v, init = n, solver = CG, eps = epsu)
            = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))
	+int2d(Th)((dx(ux)+dy(uy))*v)
	+on(1, u=sin(pi*lambda*x))
	+on(2, u=sin(pi*lambda*y))
	+on(3, u=sin(pi*lambda*x))
	+on(4, u=sin(pi*lambda*y))
	;

plot(u);

}
\end{verbatim}

2

I all a a C1 finite element in 1d (P3 hermite) it is trivial and after all is trivial.

3

a first test with the new P3 hermite finite element :

load "msh3"
load "Element_P3" // for P3
load "Element_P4" // for P3
// u_t + nu u_xxxx + u_xx + u*u_x =0; 
real L= 2*pi;
meshL Th=segment(100,[x*L]);
fespace Vh(Th,P3HL,periodic=[[1],[2]]);
fespace Ph(Th,P3);
Vh [u,ux],[v,vx],[up,upx];//  function and derivative
Ph X=x; // for 1d graphic
real nu=1; 
real dt = 0.1; 
[u,ux] = [2+sin(x),cos(x)];//  initial guess : function and derivative (need for P3 Hermite interpolation)
for (int i=0; i<100; ++i)
{
	up[]=u[]; 
     solve bilap([u,ux],[v,vx]) =
       int1d(Th)(u*v/dt +  nu*dxx(u)*dxx(v) - dx(u)*dx(v) + up*dx(u)*v) 
   - int1d(Th)(up*v/dt);
;
Ph U=u; // to plot the solution 
plot([X[],U[]], cmm=i);
}

To use this script you need to install add new P3 hermite FE in 1d,and P4 Lagrange FE in 1d .. · FreeFem/FreeFem-sources@a811767 · GitHub