Solving an advection-diffusion-like equation in a non-rectangular domain

Question

I need some guidance on solving this equation computationally for $\theta$

$$u \theta_x+v \theta_y-\nabla^2\theta=u$$

where $u$ and $v$ are known velocity functions (in $x$ and $y$ respectively) for a incompressible fluid in my domain (no slip and not penetrating into the boundary) and $f_x$ and $f_y$ are derivatives in $x$ and $y$ respectively.

My domain is a perturbed channel with periodic conditions at $x=0$ and $x=1$ and a no flux condition at the flat wall at $y=0$ but a curved wall at $y=h(x)$ s.t. $h(x)=1-a \cos(2\pi x)$ ($a$ is small let's say around 0.1).

Transforming the domain via $Y=y/H(x)$ gives me a rectangular domain, but produces 10-15 terms (which is fine but now the spacing in not uniform in the original coordinates) but at least the no flux condition is "nice" looking on the boundary.

Leaving it the same and then approximating the domain also works, but the boundary condition is a little uglier.

In either case, I am not sure what type of error this accumulates and what method to use to solve. Should I use a finite difference method in matlab? What would you recommend?

Here is a photo of the velocity field shifted for reference:

Answer 1

So here is my best attempt at a solution using Freefem, it still needs work (the mean function isn't working) and needs to be validated, but I think it can be helpful. As soon as I validate it and/or fix it, I will edit this response and remove this comment. If you know what is wrong with it, please comment:

// Parameters
real a = 0.1;   //amplitude of oscillation
real pe = 10.0; //Peclet number
int Wall = 1;
int Per1 = 2;
int Per2 = 3;
int Res = 40;   //Resolution of Solver

// Define mesh boundaries
border Gamma1(t=1, 0){x=t; y=1-a*cos(2*pi*t);label=Wall;}
border Gamma2(t=1-a, 0){x=0; y=t;label=Per1;}
border Gamma3(t=0, 1){x=t; y=0;label=Wall;}
border Gamma4(t=0, 1-a){x=1; y=t;label=Per2;}


// The triangulated domain Th is on the left side of its boundary
mesh Th = buildmesh(Gamma1(Res)+Gamma2(Res)+Gamma3(Res)+Gamma4(Res));

// Plot the Mesh
plot(Th, wait=true);

// x,y comp of velocity field
func gun0=1+a*cos(2*pi*x);
func gun1=-2*pi*a*sin(2*pi*x);
func BB=y/gun0;
func Ux=pe*(6/gun0*(BB-BB^2));
func Uy=pe*(6*gun1/gun0*BB*(BB-BB^2));

// The finite element space defined over Th is called here Vh
// Added periodic conditions on the left and right side
fespace Vh(Th, P1, periodic=[[Per1, y],[Per2, y]]);
Vh theta, v;// Define theta and v as piecewise-P1 continuous functions

//Weak form separated into bilinear and linear parts
solve Poisson(theta, v, solver=LU)
= int2d(Th)(dx(theta)*dx(v)+ dy(theta)*dy(v)
-Ux*theta*dx(v)-Uy*theta*dy(v)) // The bilinear part
-int2d(Th)(Ux*v); // The linear part

//The default boundary condition is the Neumann one (this is a direct consequence of the weak formulation)
//see https://math.stackexchange.com/questions/1537588/freefem-code-approximating-the-laplace-equation

// Output ------------------------------------------------------------
//plot of theta, the average of theta over the domain
plot(theta);
real Avg=mean(theta);

cout << " Mean=" << Avg << endl;

// Get the clock in second
real cpu=clock();

// Display the total computational time
cout << "CPU time = " << (clock()-cpu) << endl;