I wrote code on MATLAB to solve Stokes equations using Bramble and Pasciak method but I got strange values of error for pressure p

%Define problem parameters
mu = 1.0;     % viscosity
%f = [1; 0; 1; 2];   % forcing term for velocity
%g = [2; 1];        % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);

Unrecognized function or variable 'squaremesh'.

subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);
%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000; 
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time 
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H 
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B 
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B';  % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1);  % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1))  2*(eqn.A)\eqn.B_transpose; eqn.B  2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A  zeros(size(transpose(eqn.B))); zeros(size(eqn.B))  eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
    if H1ErrUCg(i) == 0 && i == 1
        H1ErrUCg(i) = max(H1ErrUCg);
    elseif H1ErrUCg(i) == 0
        H1ErrUCg(i) = H1ErrUCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrPCg(i) == 0 && i == 1
        L2ErrPCg(i) = max(L2ErrPCg);
    elseif L2ErrPCg(i) == 0
        L2ErrPCg(i) = L2ErrPCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrGradUandPCg(i) == 0 && i == 1
        L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
    elseif L2ErrGradUandPCg(i) == 0
        L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
    end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;
function [inner_pro] = inner_x_y(C, x,y)
%C = [A 0 ; 0 1]
%(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
%trans(y_1)*x_2
inner_pro = transpose(y) * C * x;
end

[Walter Roberson]

squaremesh() is missing, so we cannot test your code.

[Mostafa Kassoum]

The sequare mesh comes from the following link;

https://lyc102.github.io/ifem/projects/fem/

[Mostafa Kassoum]

Also to test the code you need the following code:

function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
z_old = x;
p_old = F - M*z_old;
r_old=p_old;
for i = 1:maxiter
    alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
    z_new   = z_old + alpha_i*p_old;
    r_new   = F - M*z_new;
    norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
    H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
    L2ErrP(i) = norm(p - z_new(length(u)+1:end));
    L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
    %L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
    relres(i) = norm_r;
    fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
        length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
    if (norm_r)< tol
        break
    end
    beta_i   = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
    p_new    = r_new - beta_i*p_old;
    % Update
    z_old = z_new;
    r_old = r_new;
    p_old = p_new;
end
end

[Walter Roberson]

The given archive does not include ExtractBoundaryEdges or StokesModelMustafaNew

%bring in external tools
sqm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/squaremesh.m";
shm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m";
unir_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/uniformrefine.m";
myu_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m";

tn = tempname();
sqm_localname = fullfile(tn, "squaremesh.m");
shm_localname = fullfile(tn, "showmesh.m");
unir_localname = fullfile(tn, "uniformrefine.m");
myu_localname = fullfile(tn, "myunique.m");

mkdir(tn);
addpath(tn)

websave(sqm_localname, sqm_url);
websave(shm_localname, shm_url);
websave(unir_localname, unir_url);
websave(myu_localname, myu_url);
%done bringing in external tools


%Define problem parameters
mu = 1.0;     % viscosity
%f = [1; 0; 1; 2];   % forcing term for velocity
%g = [2; 1];        % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);

%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);

Unrecognized function or variable 'ExtractBoundaryEdges'.

%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000; 
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time 
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H 
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B 
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B';  % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1);  % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1))  2*(eqn.A)\eqn.B_transpose; eqn.B  2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A  zeros(size(transpose(eqn.B))); zeros(size(eqn.B))  eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
    if H1ErrUCg(i) == 0 && i == 1
        H1ErrUCg(i) = max(H1ErrUCg);
    elseif H1ErrUCg(i) == 0
        H1ErrUCg(i) = H1ErrUCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrPCg(i) == 0 && i == 1
        L2ErrPCg(i) = max(L2ErrPCg);
    elseif L2ErrPCg(i) == 0
        L2ErrPCg(i) = L2ErrPCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrGradUandPCg(i) == 0 && i == 1
        L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
    elseif L2ErrGradUandPCg(i) == 0
        L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
    end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;

function [inner_pro] = inner_x_y(C, x,y)
    %C = [A 0 ; 0 1]
    %(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
    %trans(y_1)*x_2
    inner_pro = transpose(y) * C * x;
end

function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
    z_old = x;
    p_old = F - M*z_old;
    r_old=p_old;
    for i = 1:maxiter
        alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
        z_new   = z_old + alpha_i*p_old;
        r_new   = F - M*z_new;
        norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
        H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
        L2ErrP(i) = norm(p - z_new(length(u)+1:end));
        L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
        %L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
        relres(i) = norm_r;
        fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
            length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
        if (norm_r)< tol
            break
        end
        beta_i   = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
        p_new    = r_new - beta_i*p_old;
        % Update
        z_old = z_new;
        r_old = r_new;
        p_old = p_new;
    end
end

[Mostafa Kassoum]

The following code is the StokesModelMustafaNew

function pde = StokesModelMustafaNew
%% Custom Stokes equation model 
% 
pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);

% exact velocity
    function z = exactu(p)
        x = p(:,1); 
        y = p(:,2);
        z(:,1) = 1/(5*pi()^2).*sin(pi()*x).*sin(2*pi()*y);
        z(:,2) = 1/(5*pi()^2).*sin(2*pi()*x).*sin(pi()*y);
    end
% exact pressure
    function z = exactp(p)
        x = p(:,1); y = p(:,2);
        z = 2/3-x.^2-y.^2;
    end
% Dirichlet boundary condition of velocity
    function z = g_D(p)
        z = exactu(p);
    end
    function z = g_N(p)
        z(:,1) = 2*ones(size(p,1),1);
        z(:,2) = zeros(size(p,1),1);
    end
end

[Mostafa Kassoum]

I think now you can run the code and see why is the L^2 error of pressure p big value.

[Walter Roberson]

Still need ExtractBoundaryEdges

%bring in external tools
sqm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/squaremesh.m";
shm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m";
unir_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/uniformrefine.m";
myu_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m";

tn = tempname();
sqm_localname = fullfile(tn, "squaremesh.m");
shm_localname = fullfile(tn, "showmesh.m");
unir_localname = fullfile(tn, "uniformrefine.m");
myu_localname = fullfile(tn, "myunique.m");

mkdir(tn);
addpath(tn)

websave(sqm_localname, sqm_url);
websave(shm_localname, shm_url);
websave(unir_localname, unir_url);
websave(myu_localname, myu_url);
%done bringing in external tools


%Define problem parameters
mu = 1.0;     % viscosity
%f = [1; 0; 1; 2];   % forcing term for velocity
%g = [2; 1];        % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);

%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);

Unrecognized function or variable 'ExtractBoundaryEdges'.

%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000; 
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time 
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);
%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H 
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B 
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B';  % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1);  % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1))  2*(eqn.A)\eqn.B_transpose; eqn.B  2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A  zeros(size(transpose(eqn.B))); zeros(size(eqn.B))  eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
    if H1ErrUCg(i) == 0 && i == 1
        H1ErrUCg(i) = max(H1ErrUCg);
    elseif H1ErrUCg(i) == 0
        H1ErrUCg(i) = H1ErrUCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrPCg(i) == 0 && i == 1
        L2ErrPCg(i) = max(L2ErrPCg);
    elseif L2ErrPCg(i) == 0
        L2ErrPCg(i) = L2ErrPCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrGradUandPCg(i) == 0 && i == 1
        L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
    elseif L2ErrGradUandPCg(i) == 0
        L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
    end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;

function [inner_pro] = inner_x_y(C, x,y)
    %C = [A 0 ; 0 1]
    %(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
    %trans(y_1)*x_2
    inner_pro = transpose(y) * C * x;
end

function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
    z_old = x;
    p_old = F - M*z_old;
    r_old=p_old;
    for i = 1:maxiter
        alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
        z_new   = z_old + alpha_i*p_old;
        r_new   = F - M*z_new;
        norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
        H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
        L2ErrP(i) = norm(p - z_new(length(u)+1:end));
        L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
        %L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
        relres(i) = norm_r;
        fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
            length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
        if (norm_r)< tol
            break
        end
        beta_i   = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
        p_new    = r_new - beta_i*p_old;
        % Update
        z_old = z_new;
        r_old = r_new;
        p_old = p_new;
    end
end

function pde = StokesModelMustafaNew
    %% Custom Stokes equation model 
    % 
    pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);
    
    % exact velocity
        function z = exactu(p)
            x = p(:,1); 
            y = p(:,2);
            z(:,1) = 1/(5*pi()^2).*sin(pi()*x).*sin(2*pi()*y);
            z(:,2) = 1/(5*pi()^2).*sin(2*pi()*x).*sin(pi()*y);
        end
    % exact pressure
        function z = exactp(p)
            x = p(:,1); y = p(:,2);
            z = 2/3-x.^2-y.^2;
        end
    % Dirichlet boundary condition of velocity
        function z = g_D(p)
            z = exactu(p);
        end
        function z = g_N(p)
            z(:,1) = 2*ones(size(p,1),1);
            z(:,2) = zeros(size(p,1),1);
        end
end

[Mostafa Kassoum]

The following code is ExtractBoundaryEdges.

function [Dirichlet,bdFlag] = ExtractBoundaryEdges(elem, node)

bdFlag = setboundary(node,elem,'Dirichlet','(x==-1) | (x==1) | (y==1) | (y==-1)');%,'Neumann', '(y==1)');

% bdFlag = setboundary(node,elem,'Neumann','(y==1)');

% [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);

% [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);

allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];

Dirichlet = allEdge((bdFlag(:) == 1),:);

end

[Walter Roberson]

Need StokesisoP2P0

%bring in external tools
sqm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/squaremesh.m";
shm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m";
unir_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/uniformrefine.m";
myu_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m";
setb_url = "https://raw.githubusercontent.com/lyc102/ifem/master/fem/setboundary.m";

tn = tempname();
sqm_localname = fullfile(tn, "squaremesh.m");
shm_localname = fullfile(tn, "showmesh.m");
unir_localname = fullfile(tn, "uniformrefine.m");
myu_localname = fullfile(tn, "myunique.m");
setb_localname = fullfile(tn, "setboundary.m");

mkdir(tn);
addpath(tn)

websave(sqm_localname, sqm_url);
websave(shm_localname, shm_url);
websave(unir_localname, unir_url);
websave(myu_localname, myu_url);
websave(setb_localname, setb_url);
%done bringing in external tools


%Define problem parameters
mu = 1.0;     % viscosity
%f = [1; 0; 1; 2];   % forcing term for velocity
%g = [2; 1];        % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);

%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000; 
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time 
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);

Unrecognized function or variable 'StokesisoP2P0'.

%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H 
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B 
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B';  % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1);  % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1))  2*(eqn.A)\eqn.B_transpose; eqn.B  2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A  zeros(size(transpose(eqn.B))); zeros(size(eqn.B))  eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
    if H1ErrUCg(i) == 0 && i == 1
        H1ErrUCg(i) = max(H1ErrUCg);
    elseif H1ErrUCg(i) == 0
        H1ErrUCg(i) = H1ErrUCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrPCg(i) == 0 && i == 1
        L2ErrPCg(i) = max(L2ErrPCg);
    elseif L2ErrPCg(i) == 0
        L2ErrPCg(i) = L2ErrPCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrGradUandPCg(i) == 0 && i == 1
        L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
    elseif L2ErrGradUandPCg(i) == 0
        L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
    end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;

function [inner_pro] = inner_x_y(C, x,y)
    %C = [A 0 ; 0 1]
    %(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
    %trans(y_1)*x_2
    inner_pro = transpose(y) * C * x;
end

function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
    z_old = x;
    p_old = F - M*z_old;
    r_old=p_old;
    for i = 1:maxiter
        alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
        z_new   = z_old + alpha_i*p_old;
        r_new   = F - M*z_new;
        norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
        H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
        L2ErrP(i) = norm(p - z_new(length(u)+1:end));
        L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
        %L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
        relres(i) = norm_r;
        fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
            length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
        if (norm_r)< tol
            break
        end
        beta_i   = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
        p_new    = r_new - beta_i*p_old;
        % Update
        z_old = z_new;
        r_old = r_new;
        p_old = p_new;
    end
end

function pde = StokesModelMustafaNew
    %% Custom Stokes equation model 
    % 
    pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);
    
    % exact velocity
        function z = exactu(p)
            x = p(:,1); 
            y = p(:,2);
            z(:,1) = 1/(5*pi()^2).*sin(pi()*x).*sin(2*pi()*y);
            z(:,2) = 1/(5*pi()^2).*sin(2*pi()*x).*sin(pi()*y);
        end
    % exact pressure
        function z = exactp(p)
            x = p(:,1); y = p(:,2);
            z = 2/3-x.^2-y.^2;
        end
    % Dirichlet boundary condition of velocity
        function z = g_D(p)
            z = exactu(p);
        end
        function z = g_N(p)
            z(:,1) = 2*ones(size(p,1),1);
            z(:,2) = zeros(size(p,1),1);
        end
end

function [Dirichlet,bdFlag] = ExtractBoundaryEdges(elem, node)
        bdFlag = setboundary(node,elem,'Dirichlet','(x==-1)  | (x==1) | (y==1) | (y==-1)');%,'Neumann', '(y==1)');
    %     bdFlag = setboundary(node,elem,'Neumann','(y==1)');
        %     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
    %     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
        allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];
        Dirichlet = allEdge((bdFlag(:) == 1),:);
    
end

[Mostafa Kassoum]

You need download the codes in the following link on MATLAB to get what you need.

https://lyc102.github.io/ifem/projects/fem/

[Walter Roberson]

%bring in external tools
sqm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/squaremesh.m";
shm_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/showmesh.m";
unir_url = "https://raw.githubusercontent.com/lyc102/ifem/master/mesh/uniformrefine.m";
myu_url = "https://raw.githubusercontent.com/lyc102/ifem/master/tool/myunique.m";
setb_url = "https://raw.githubusercontent.com/lyc102/ifem/master/fem/setboundary.m";
S2P0_url = "https://raw.githubusercontent.com/lyc102/ifem/master/equation/StokesisoP2P0.m";
dof_url = "https://raw.githubusercontent.com/lyc102/ifem/master/dof/dof3P2.m";

tn = tempname();
sqm_localname = fullfile(tn, "squaremesh.m");
shm_localname = fullfile(tn, "showmesh.m");
unir_localname = fullfile(tn, "uniformrefine.m");
myu_localname = fullfile(tn, "myunique.m");
setb_localname = fullfile(tn, "setboundary.m");
S2P0_localname = fullfile(tn, "StokesisoP2P0.m");
dof3_localname = fullfile(tn, "dof3P2.m");
dof_localname = fullfile(tn, "dofP2.m");

mkdir(tn);
addpath(tn)

websave(sqm_localname, sqm_url);
websave(shm_localname, shm_url);
websave(unir_localname, unir_url);
websave(myu_localname, myu_url);
websave(setb_localname, setb_url);
websave(S2P0_localname, S2P0_url);
websave(dof_localname, dof_url);
websave(dof3_localname, dof_url);

%done bringing in external tools


%Define problem parameters
mu = 1.0;     % viscosity
%f = [1; 0; 1; 2];   % forcing term for velocity
%g = [2; 1];        % forcing term for pressure
%% Mesh creation
h = 0.5; % discretization step
[node,elem] = squaremesh([-1,1,-1,1],h);
subplot(1,2,1); showmesh(node,elem);
N = size(node,1);
NTc = size(elem,1);
% Uniform refinement to get a fine grid
[node,elem] = uniformrefine(node,elem);
% Uniform refinement to get a fine grid
subplot(1,2,2); showmesh(node,elem);

%%extract boundary conditions
[Dirichlet, bdFlag] = ExtractBoundaryEdges(elem, node);
%% Load Stokes Data Partial Differential Equations
%Different velocity and pressure equation models
pde = StokesModelMustafaNew;%StokesModelMustafa;%StokesModelMustafaNew;%Stokesdata3;%StokesModelMustafa2;%StokesModelMustafaSin;%;
option.elemType = 'isoP2P0';
option.fquadorder = 3;
tolerance = 1e-3;
max_iterations = 1000; 
mesh.node = node;
mesh.elem = elem;
mesh.bdFlag = bdFlag;
%The following line of code is only run to get the eqn (equation structure) easier without defining it every time 
[soln,eqn,info] = StokesisoP2P0(node,elem,bdFlag,pde, option);

Index in position 2 exceeds array bounds. Index must not exceed 3.

Error in dofP2 (line 17)
totalEdge = uint32([elem(:,[1 2]); elem(:,[1 3]); elem(:,[1 4]); ...

Error in StokesisoP2P0 (line 25)
[tempvar,edgeC] = dofP2(elemC);

%% Call the solvers and calulate L2 errors, run time and outcome as well as convergence rate
% [soln,eqn,err,time,solver] = femStokes(mesh,pde,option);
% return
uI = pde.exactu([node; (node(eqn.edge(:,1),:)+node(eqn.edge(:,2),:))/2]);
uI = reshape(uI, [2*length(uI),1]);
pI = pde.exactp([node(elem(:,1),1) node(elem(:,3),2)]);
% Define matrix H 
%H = [1 0 0 2; 0 2 1 0; 2 3 1 0; 0 1 0 1];
% Define matrix A
%A = H * H';
% Define matrix B 
%B = [1 2 1 0; 0 1 0 0];
% Assemble the finite element system
%[A, B, ~, ~] = assempde(pde,node,elem,'f',f,'g',g,'mu',mu);
tic;
% Assemble the Stokes matrix
N = size(eqn.A,1);
M = size(eqn.B,1);
%StokesMatrix = [eqn.A, eqn.B'; eqn.B, sparse(M,M)];
% Define dimensions of eqn.B'
eqn.B_transpose = eqn.B';  % Transpose of eqn.B
% Construct a sparse MxM matrix (adjust M to your desired size)
M = size(eqn.B, 1);  % Number of columns in eqn.B (or rows in eqn.B')
sparse_M_M = sparse(M, M);
% Check dimensions of matrices involved in the problematic line
size_eqn.B = size(eqn.B);
size_eqn.A = size(eqn.A);
% Perform the matrix multiplication explicitly and verify dimensions
result = 2 * eqn.B / (eqn.A)* eqn.B_transpose;
size_result = size(result);
% Concatenate eqn.A, eqn.B_transpose, eqn.B, and sparse_M_M
StokesMatrix = [eqn.A, eqn.B_transpose; eqn.B, sparse_M_M];
% Define matrix M
%M = [2*eye(size(eqn.A,1))  2*(eqn.A)\eqn.B_transpose; eqn.B  2*eqn.B*(eqn.A)\eqn.B_transpose];
M_top = [2 * eye(size(eqn.A, 1)), 2 * (eqn.A\eqn.B_transpose)];
M_bottom = [eqn.B, 2 * eqn.B * (eqn.A \ eqn.B_transpose)];
M = [M_top; M_bottom];
% Define matrix F
F = [2 * eqn.A \ eqn.f; 2 * eqn.B * (eqn.A \ eqn.f) - eqn.g];
% Define the right-hand side
rhs = [eqn.f; eqn.g];
x0 = ones(size(F));
% Solve the Stokes system
solution = StokesMatrix \ rhs;
sol = M \ F;
% Extract the velocity and pressure components from the solution
u = solution(1:N);
p = solution(N+1:end);
% Calculate the residual R
R = F - M * sol;
% Define matrix C
C = [0.5 * eqn.A  zeros(size(transpose(eqn.B))); zeros(size(eqn.B))  eye(size(eqn.B,1))];
% Define an appropriate non-zero vector x
[res_conjgrad, iterCg, relresCg, H1ErrUCg, L2ErrPCg, L2ErrGradUandPCg] = conjgradNew3(eqn.A, C,F ,[eqn.f; eqn.g],x0 , max_iterations, tolerance, u, p, h, M);
for i = 1:iterCg
    if H1ErrUCg(i) == 0 && i == 1
        H1ErrUCg(i) = max(H1ErrUCg);
    elseif H1ErrUCg(i) == 0
        H1ErrUCg(i) = H1ErrUCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrPCg(i) == 0 && i == 1
        L2ErrPCg(i) = max(L2ErrPCg);
    elseif L2ErrPCg(i) == 0
        L2ErrPCg(i) = L2ErrPCg(i-1);
    end
end
for i = 1:iterCg
    if L2ErrGradUandPCg(i) == 0 && i == 1
        L2ErrGradUandPCg(i) = max(L2ErrGradUandPCg);
    elseif L2ErrGradUandPCg(i) == 0
        L2ErrGradUandPCg(i) = L2ErrGradUandPCg(i-1);
    end
end
p = res_conjgrad(size(eqn.A,1)+1:end);
u = res_conjgrad(1:size(eqn.A,1));
time_BPCG_finish = toc;

function [inner_pro] = inner_x_y(C, x,y)
    %C = [A 0 ; 0 1]
    %(x,y) = (Cx,y) = (Ax_1,y_1) + (x_2,y_2) = transpose(x_1)*A*y_1 +transpose(x_2)*y_2
    %trans(y_1)*x_2
    inner_pro = transpose(y) * C * x;
end

function [x, i, relres, H1ErrU, L2ErrP, L2ErrGradUandP] = conjgradNew3(~, C,F,~, x, maxiter, tol, u, p, h, M)
    z_old = x;
    p_old = F - M*z_old;
    r_old=p_old;
    for i = 1:maxiter
        alpha_i = inner_x_y(C, r_old,p_old) / inner_x_y(C, M*p_old,p_old);
        z_new   = z_old + alpha_i*p_old;
        r_new   = F - M*z_new;
        norm_r = transpose(r_new)*r_new;%inner_x_y(C, r_new,r_new);
        H1ErrU(i) = norm(gradient(u - z_new(1:length(u))));
        L2ErrP(i) = norm(p - z_new(length(u)+1:end));
        L2ErrGradUandP(i) = sqrt((H1ErrU(i))^2 + (L2ErrP(i))^2);
        %L2ErrGradUandP(i) = sqrt((norm(gradient(z_new(1:length(uI)),h) - gradient(uI,h)))^2 + (norm(pI - z_new(length(uI)+1:end)))^2);
        relres(i) = norm_r;
        fprintf('#dof: %8.0u, Bramble Pasciak CG iter: %2.0u, err = %12.8g\n \n',...
            length(z_new)+length(z_new(length(u)+1:end)), i, norm_r);
        if (norm_r)< tol
            break
        end
        beta_i   = inner_x_y(C, M*r_new,p_old) / inner_x_y(C, M*p_old,p_old);
        p_new    = r_new - beta_i*p_old;
        % Update
        z_old = z_new;
        r_old = r_new;
        p_old = p_new;
    end
end

function pde = StokesModelMustafaNew
    %% Custom Stokes equation model 
    % 
    pde = struct('f', 0, 'exactp', @exactp, 'exactu', @exactu,'g_D',@g_D,'g_N',@g_N);
    
    % exact velocity
        function z = exactu(p)
            x = p(:,1); 
            y = p(:,2);
            z(:,1) = 1/(5*pi()^2).*sin(pi()*x).*sin(2*pi()*y);
            z(:,2) = 1/(5*pi()^2).*sin(2*pi()*x).*sin(pi()*y);
        end
    % exact pressure
        function z = exactp(p)
            x = p(:,1); y = p(:,2);
            z = 2/3-x.^2-y.^2;
        end
    % Dirichlet boundary condition of velocity
        function z = g_D(p)
            z = exactu(p);
        end
        function z = g_N(p)
            z(:,1) = 2*ones(size(p,1),1);
            z(:,2) = zeros(size(p,1),1);
        end
end

function [Dirichlet,bdFlag] = ExtractBoundaryEdges(elem, node)
        bdFlag = setboundary(node,elem,'Dirichlet','(x==-1)  | (x==1) | (y==1) | (y==-1)');%,'Neumann', '(y==1)');
    %     bdFlag = setboundary(node,elem,'Neumann','(y==1)');
        %     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
    %     [node,elem,bdFlag] = uniformbisect(node,elem,bdFlag);
        allEdge = [elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])];
        Dirichlet = allEdge((bdFlag(:) == 1),:);
    
end

[Mostafa Kassoum]

Did the code run?

[Torsten]

@Mostafa Kassoum

No. You see it from the error message above in red.

If you include code here in the forum, you can (and should) test if it's at least technically working by clicking the green Run button in the menu bar.

[Mostafa Kassoum]

It worked for me anf I got the results for everything and they are nice results except the error of pressure P, they are so big. I think there is a mistake but I didn't got it.

[Mostafa Kassoum]

Foe example I got the following results:

H1 error for velocity u are: 0.129922159092270 0.138027304457401

L2 error for pressure p are: 12.6454430756918 12.6450517755845

I need know why are L2 error p big

[Walter Roberson]

You can copy my entire post to your system and test it locally.

I suspect that something you have installed does not match the latest code at https://github.com/lyc102/ifem

[Mostafa Kassoum]

I copied your entire post and runned it but I got also a big value of L2 error 11.3255614294248 11.3251245255107

As you can see they are almost equal to the values I obtained previously. Where do you think the mistake is?