[fedvasu]

Hi Foamers,

I am using a simple internal faces only scheme for smoothing, I can do this just fine in serial.

Basically I am by hand setting coefficients of the matrices.

my system is Gx=W, I have verified that this is exactly what I need.

    fvScalarMatrix W(fvm::Sp(rho,rhoNew)); //my rhs matrix
    fvScalarMatrix G(fvm::Sp(rho,rhoNew)); //lhs matrix


    W.diag() = rho.primitiveField();

    const labelUList& owner = mesh.owner(); //internal faces 
    const labelUList& neighbour = mesh.neighbour(); // internal faces neighbour

    forAll(owner,facei)
    {
        G.upper()[facei] += (-alpha*(rho[owner[facei]] + rho[neighbour[facei]])/2); // r[facei] = (r[owner] + r[neighbour]) /2;
    }

    G.negSumDiag();
    G.source() += W.diag()*rho.primitiveField(); //rhs 

    G.relax();
    G.solve();

I would like to know how to handle special procBoundary "boundary patches", how should I index my G and W

[fedvasu]

Originally Posted by fedvasu

Hi Foamers,

I am using a simple internal faces only scheme for smoothing, I can do this just fine in serial.

Basically I am by hand setting coefficients of the matrices.

my system is Gx=W, I have verified that this is exactly what I need.

    fvScalarMatrix W(fvm::Sp(rho,rhoNew)); //my rhs matrix
    fvScalarMatrix G(fvm::Sp(rho,rhoNew)); //lhs matrix


    W.diag() = rho.primitiveField();

    const labelUList& owner = mesh.owner(); //internal faces 
    const labelUList& neighbour = mesh.neighbour(); // internal faces neighbour

    forAll(owner,facei)
    {
        G.upper()[facei] += (-alpha*(rho[owner[facei]] + rho[neighbour[facei]])/2); // r[facei] = (r[owner] + r[neighbour]) /2;
    }

    G.negSumDiag();
    G.source() += W.diag()*rho.primitiveField(); //rhs 

    G.relax();
    G.solve();

I would like to know how to handle special procBoundary "boundary patches", how should I index my G and W

This took 4-5 days for me to actually figure out, but this book is quite helpful
https://www.springer.com/us/book/9783319168739

    fvScalarMatrix W(fvm::Sp(rho,rhoNew)); //my rhs matrix
    fvScalarMatrix G(fvm::Sp(rho,rhoNew)); //lhs matrix


    W.diag() = rho.primitiveField();

    const labelUList& owner = mesh.owner(); //internal faces 
    const labelUList& neighbour = mesh.neighbour(); // internal faces neighbour

   //this instead of doing the main loop by hand
   surfaceScalarField rhoS("rhoS",fvc::interpolate(rho,"midPoint"));
    A.upper() = -alphaS*rhoS;
    A.negSumDiag();
  


    G.diag() = A.diag() + W.diag();
    G.upper() = A.upper();

    //works
    G.source() = W.diag()*sfield.primitiveField();

    //this is how you handle parallelism using boundaryCoeffs and internallCoeffs on processor patches
    forAll(mesh_.boundary(), patchi)
    {
        const polyPatch &ppatch = mesh_.boundaryMesh()[patchi];

        if (isA<processorPolyPatch>(ppatch))
        {
            auto rhoSBf = rhoS.boundaryFieldRef()[patchi];
            G.internalCoeffs()[patchi] = alphaS*rhoSBf;
            G.boundaryCoeffs()[patchi] = alphaS*rhoSBf;
        }
    }

    G.relax();
    G.solve();

This specific issue is solved for me.

[dlahaye]

Thanks for sharing.
Domenico Lahaye.