Smoluchowski / Diffusion Drift Equation Help and Adaptive Meshing

General Discussion
1

Dear All,

I am trying to solve the diffusion equation between two parallel plates in the x direction. This problem is complicated by the fact that these diffusive particles are “active”, so they are orientable and have a drift force in the direction of their orientation (which is the y variable in this code).

An issue that I am having is that it seems it is very sensitive to mesh size, even when there should not be any sharp spikes like one would see in a boundary layer. Do you have any advice on how to make equations more robust to mesh size, and is this a problem common to the Diffusion/Drift Equation?

Lastly, because I am using no flux and periodic boundary conditions, I need to use Lagrange multipliers to choose one particular solution, would that choice affect the sensitivity of the meshing, and how does one implement adaptive meshing for something with varf and lagrange multipliers? I do not see a clear example in the documentation. https://doc.freefem.org/examples/mesh-generation.html#mesh-adaptation

Thank you, I have attached a basic version of my code below.

load "MUMPS_seq"
searchMethod=1; // More safe seach algo

////// Parameters
/// Lagrange Multiplier
real Nintegration = 1.0;

//INSERTVARIABLES
//START VARIABLE DEFS
real xBoxEdgeLen=5.0;
real difft=1.0;
real taur=1.0;
real u0=1.0;
real Lc=1.0;
// END VARIABLE DEFS 



int nnxlen = 100, ny = 100;

real h = 2*pi;
func ymin = 0;
func ymax = h;
real ax = xBoxEdgeLen;
// Mesh 2D

mesh Th = square(nnxlen, ny, [-ax/2+(ax)*x, (ymax)*y]);;
{ // Cleaning the memory correctly
    int[int] old2new(0:Th.nv-1);
    fespace Vh2(Th, P1);
    Vh2 sorder = x + y;
    sort(sorder[], old2new);
    int[int] new2old = old2new^-1; // Inverse permutation
    Th = change(Th, renumv=new2old);
    sorder[] = 0:Th.nv-1;
}
{
    fespace Vh2(Th, P1);
    Vh2 nu;
    nu[] = 0:Th.nv-1;
    //plot(nu, cmm="nu=", wait=true);
}


real pes = u0*Lc/difft, gammasqrd = (Lc*Lc)/(difft*taur);
fespace Vh(Th, P1, periodic=[[1, x], [3, x]]);
int n = Vh.ndof;
int n1 = n+1;
Vh probability, testfunc; // probabiility and test function.

// Problem
varf va(probability,testfunc) = // definion of the problem
  int2d(Th)(dx(testfunc)*dx(probability) +  gammasqrd*dy(probability)*dy(testfunc))
  +int2d(Th)(testfunc*pes*(cos(y)*dx(probability)) )
  +int1d(Th,2)(1*pes*(probability*testfunc*(cos(y))))
  +int1d(Th,4)(-1*pes*(probability*testfunc*(cos(y))))
;

func f = 0;//
varf vL (probability, testfunc) = int2d(Th)(f*testfunc);
varf vb (probability, testfunc) = int2d(Th)(1.*testfunc); /// Multiplying this by 10 is like *dividing* b1 by 10.

matrix A = va(Vh, Vh);
real[int] b = vL(0, Vh);
real[int] B = vb(0, Vh);

// Block matrix
matrix AA = [ [ A, B ], [ B', 0 ] ];
set(AA, solver=sparsesolver);

real[int] bb(n+1), xx(n+1), b1(1), l(1);
b1 = Nintegration; // So I think this is saying what the function integrates too. I.e. int2d(probability) == b1?
// Builds the right hand side block
bb = [b, b1];
//bb =[0,b1];
// Solve
xx = AA^-1 * bb;

// Set values
[probability[],l] = xx;

// Display
cout << " l = " << l(0) << " , b(u, 1) =" << B'*probability[] << endl;

// Plot
//plot(Th,wait=true);
plot(probability);


// // Here I just move the solutions to a data file so I can analyze it.
ofstream file("TwoWalls_SL_trial.dat");
for (int i = 0; i < Th.nv; i++) {
    real xCoord = Th(i).x;
    real yCoord = Th(i).y;

    real val = probability(xCoord, yCoord);
    file << xCoord << " " << yCoord << " " << val << endl;
}