[mahdi_kh8]

I am trying to check different pressure solvers and have better understanding.

Interestingly, when I replace codes I do not get correct results.

Original code for pressure:

    Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye(nx));
    Lp(1,1) = 3/2*Lp(1,1);
    perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp';
    rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1);
    p(perp) = -Rp\(Rpt\rhs(perp));
    P = reshape(p,nx,ny);

Edited code:

    Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye(nx));
    Lp(1,1) = 3/2*Lp(1,1);
    perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp';
    rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1);
    % p(perp) = -Rp\(Rpt\rhs(perp));

    p= Lp\rhs;

    P = reshape(p,nx,ny);

or

    Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye(nx));
    Lp(1,1) = 3/2*Lp(1,1);
    perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp';
    rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1);
    % p(perp) = -Rp\(Rpt\rhs(perp));

    [p,flagCG,relresCG,iterCG,rCG] = pcg(Lp,rhs,1e-06,10000);

    P = reshape(p,nx,ny);

None of the above changes work.

What is missing here?

Any advise would be appreciated.

Cheers,

[mahdi_kh8]

There was a small mistake:

    Lp = kron(speye(ny),K1(nx,hx,1))+kron(K1(ny,hy,1),speye(nx));
    Lp(1,1) = 3/2*Lp(1,1);
    perp = symamd(Lp); Rp = chol(Lp(perp,perp)); Rpt = Rp';
    rhs = reshape(diff([uW;U;uE])/hx+diff([vS' V vN']')'/hy,[],1);
    % p(perp) = -Rp\(Rpt\rhs(perp));

    [p,flagCG,relresCG,iterCG,rCG] = pcg(Lp,rhs,1e-06,10000);

    P = reshape(p,nx,ny);
    P=-P;