Currently struggling with boundary conditions for 2d advection equation.

left dirichlet u = (1,0)
bot top dirichlet u = (0,0)
right neumann du/dx = (0,0)

topf botf leftf rightf - arrays that contain indicies of boundary faces
nbfaces - # of boundary faces
uface - array of vectors containing u and v at faces (inner ones are interpolated)
ucell - array of vectors containing u and v at cells
unode - array of vectors containing u and v at nodes
nnodes - nu of nodes
link_... - array for link data

% interpolate inner nodes
unode = zeros(nnodes,2);
for iv = 1:nnodes
    unode(iv,:) = 0;
    for ic = 1:sum(wv(iv,:)~=0) % over nonzero wv elements
        if bnode(iv) == 0 % for interior nodes
            unode(iv,:) = unode(iv,:) + ucell(link_node_to_cell(iv,ic),:) * wv(iv,ic);
% hardwire left boundary nodes
for ifc = 1:nbfaces
    if ismember(ifc,leftf)
        unode(link_face_to_node(ifc,1),1) = 1;
        unode(link_face_to_node(ifc,2),1) = 1;

% interpolation over inner faces 
for ifc = 1:nfaces
    c1 = link_face_to_cell(ifc,1);
    c2 = link_face_to_cell(ifc,2);
    if c2 ~= 0
        uface(ifc,:) = wf(ifc) * ucell(c1,:) + (1-wf(ifc)) * ucell(c2,:);

managed to to the most silly thing - nullify the cells that have bottom boundary faces

% hardwire boundary faces
for ifc = 1:nbfaces
    if ismember(ifc,botf)
        uface(ifc,1) = 0;
        ic = link_face_to_cell(ifc,1);
        ucell(ic,:) = 0;
    elseif ismember(ifc,topf)
        uface(ifc,1) = 0;
    elseif ismember(ifc,leftf)
        uface(ifc,1) = 1;
    elseif ismember(ifc,rightf)
        ic = link_face_to_cell(ifc,1);
        uface(ifc,1) = ucell(ic,1);

advect = zeros(ncells,2); 
for ic = 1:ncells
    advect(ic,:) = 0;
    for ifc = 1:3
        fv = link_cell_to_face(ic,ifc);
        fsign = snsign(ic,ifc);
        udotn = uface(fv,:) *( sn(fv,:) * fsign)';
        gf = udotn * areaf(fv); % mass according for upwind scheme
        advect(ic,:) = advect(ic,:) + (max(gf,0) * ucell(ic,:) + min(gf,0) * uface(fv,:));

% solve the cell
diffu = zeros(ncells,2);
dt = 1e-4;
for ic = 1:ncells
    fc = link_cell_to_face(ic,1);
        ucellnext(ic,:) = dt * (diffu(ic,:) -   advect(ic,:)) / vol(ic) + ucell(ic,:);
ucell = ucellnext

Receiving weird results, that:
First of all diverge.
And my method of addressing the bottom boundary looks like crutches.

When trying the code without nullifying the boundary cells received diverging results at the bottom boundary.


Only boundary cells behaved strangely.

enter image description here

Could one recommend the literature how to address the boundaries? Or point out the reason for the issue.
Self search lead to articles with pure theoretical descriptions.

Can paste the whole code with input file if necessary. (~600 lines)

Strange results with boundary cells being nullified, the result might have converged, however the velocities at the left were set to 1 (after searching for cells with values over 1.0, a plenty of them seem to be the cells near the boundary cells at the bottom):

diverging solution?


PDE problem statement.
$$\frac{\partial }{\partial t}(u)+\nabla\cdot(uu)=0 $$ boundaries: $$u(left,t) = 1$$ $$\frac{\partial u}{\partial x}(right,t) = 0$$ $$u(top,t) = 0$$ $$u(bot,t) = 0$$ Discretisation derivation, integrate over volume:

$$\int_{V_O} \frac{\partial }{\partial t}(u)dV + \int_{V_O}\nabla\cdot(uu)dV=0$$
applying Gauss divergence theorem to the advection term:

$$\int_{V_O} \frac{\partial }{\partial t}(u)dV + \int_{S}(uu)\cdot\hat{n}dA=0$$

using upwind difference scheme: $$G_f=(u\cdot\hat{n}A)_f$$ $$\frac{u^{n+1}-u^{n}}{\Delta t}V_O+\sum_{f=1}^{N_f,O}\biggl( max(G_f,0)u_O + min(G_f,0)u_{N(f)} \biggl)=0$$


$$u^{n+1}=u^n-\frac{\Delta t}{V_O}\sum_{f=1}^{N_f,O}\biggl( max(G_f,0)u_O + min(G_f,0)u_{N(f)} \biggl)$$

boundaries separated in 4 different groups:
left faces $u_{f,L} = 1$
right faces $u_{f,R} = u_{cell,R}$
top faces $u_{f,T} = 0$
bot faces $u_{f,B} = 0$

Nodes take the values of faces that they belong to.

Inner face and nodes are interpolated from cells values based on inverse distance.

Problem arises at the bottom boundary faces and cells.

  • 1
    $\begingroup$ Can you describe your problem mathematically? Add the PDE and describe the numerical scheme? $\endgroup$ – nicoguaro Mar 26 at 23:53
  • $\begingroup$ @nicoguaro just did, thank you $\endgroup$ – 2Napasa Mar 28 at 12:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.