[danny123]

Hello,

I would like to overwrite individual matrix values in OpenFoam. I understand in order to get the values you use simply:

vectorvalue = vector[cellI] whereas vectorvalue is a scalar and vector a volScalarField in my example. Still I get some trouble to do the opposite operation:

vector[cellI] = vectorvalue

It does not seem to update, e.g. applied on alpha1 in interFoam:

        scalar prec = 1e-10;

        volScalarField alpha1n(alpha1);

        forAll(alpha1n,cellI)
        {
            if (alpha1[cellI] <= 0.0)
            {
                alpha1n[cellI] = 0.0;
            }
            else
            {
                if ((alpha1[cellI] > 0.0) and (alpha1[cellI] < prec))
                {
                    alpha1n[cellI] = 1.0 - Foam::cos(0.5*Foam::constant::mathematical::pi*alpha1[cellI]/prec);
                }
                else
                {
                    if (alpha1[cellI] >= 1.0)
                    {
                        alpha1n[cellI] = 1.0;
                    }
                    else
                    {
                        if (alpha1[cellI] > 1.0 - prec) 
                        {
                            alpha1n[cellI] = Foam::cos(0.5*Foam::constant::mathematical::pi*(1.0-alpha1[cellI])/prec);
                        }
                        else
                        {
                            alpha1n[cellI] = alpha1[cellI];
                        }
                    }
                }
            }
        }

        Info<< "Phase-1 volume fraction (routine 1) = "
            << alpha1.weightedAverage(mesh.Vsc()).value()
            << "  Min(alpha1) = " << min(alpha1).value()
            << "  Max(alpha1) = " << max(alpha1).value()
            << "  Min(alpha1n) = " << min(alpha1n).value()
            << "  Max(alpha1n) = " << max(alpha1n).value()
            << endl;
        
        alpha1 = alpha1n;

If I run this code, the initial assignement seems to work (alpha1n is initiaized to be alpha1), but then the individual value assignment does not seem to work properly (alpha 1 should be between 0 and 1):

PIMPLE: iteration 3
smoothSolver:  Solving for alpha.epdm, Initial residual = 0.000114686, Final residual = 6.61771e+107, No Iterations 1000
Phase-1 volume fraction (routine 1) = 1.08333e+113  Min(alpha1) = 0  Max(alpha1) = 5.03793e+117  Min(alpha1n) = 0  Max(alpha1n) = 5.03793e+117

Is there anybody who could help?

Related question: How do I get values of system of equation, such as fvScalarMatrix. I understand that this is linear equation system such as A x = B, A being a matrix and x and B vector each, whereas x is the unknown to be looked for. I need the values of A and B.

Regards,

Daniel

[danny123]

The code actually works as expected limiting alpha1 between 0 and 1. However the min and max values are not printed out correctly .

For the 2nd matrix problem, I found nn.source as a way to get eg. A * p_rgh = B for B:

pEqn_s = p_rghEqn.source();

The problem now is what type pEqn_s shall be. I have initialized it to be identical to p_rgh, because this sounds logical. However I get following error:

pEqn_m1.H:46:34: Fehler: keine Übereinstimmung für  »operator=« in »pEqn_s = p_rghEqn.Foam::fvMatrix<Type>::source  [with Type = double]()«
pEqn_m1.H:46:34: Anmerkung: Kandidaten sind:
/home/dw/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricField.C:1080:6:  Anmerkung: void Foam::GeometricField<Type, PatchField,  GeoMesh>::operator=(const Foam::GeometricField<Type, PatchField,  GeoMesh>&) [with Type = double, PatchField = Foam::fvPatchField,  GeoMesh = Foam::volMesh]
/home/dw/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricField.C:1080:6:  Anmerkung:   keine bekannte Umwandlung für Argument 1 von  »Foam::Field<double>« nach »const Foam::GeometricField<double,  Foam::fvPatchField, Foam::volMesh>&«
/home/dw/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricField.C:1105:6:  Anmerkung: void Foam::GeometricField<Type, PatchField,  GeoMesh>::operator=(const Foam::tmp<Foam::GeometricField<Type,  PatchField, GeoMesh> >&) [with Type = double, PatchField =  Foam::fvPatchField, GeoMesh = Foam::volMesh]
/home/dw/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricField.C:1105:6:  Anmerkung:   keine bekannte Umwandlung für Argument 1 von  »Foam::Field<double>« nach »const  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField,  Foam::volMesh> >&«
/home/dw/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricField.C:1141:6:  Anmerkung: void Foam::GeometricField<Type, PatchField,  GeoMesh>::operator=(const Foam::dimensioned<Form>&) [with  Type = double, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]
/home/dw/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/GeometricField.C:1141:6:  Anmerkung:   keine bekannte Umwandlung für Argument 1 von  »Foam::Field<double>« nach »const  Foam::dimensioned<double>&«

If A is a matrix of n x n size and p_rgh a vector having n rows, B should have n rows too. Dimensions should be different (same as p_rghEqn). In Createfields.H dimensions are not set, but by a file in the 0 directory, which does not exist yet at compilation:

    Info<< "Reading field pEqn_s\n" << endl;
    volScalarField pEqn_s
    (
        IOobject
        (
            "pEqn_s",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

Any help?

Daniel

[danny123]

Hello again,

Well there seems nobody out there, though there is time between games in World Cup. I found some solutions of my problem meanwhile, but now I am stuck again.

The solution to the type in fvScalaMatrix seems to be the following:

        Foam::Field<double> pEqn_s = p_rghEqn.source();
        Foam::Field<double> pEqn_d = p_rghEqn.diag();
        Foam::Field<double> pEqn_u = p_rghEqn.upper();
        Foam::Field<double> pEqn_l = p_rghEqn.lower();

This compiles at least, I guess it should work too. Thus I have all values of A matrix and B vector, Now I need the address of the upper and lower triangle. I tried the following:

        Foam::Field<int> pEqn_u_addr = p_rghEqn.upperAddr();
        Foam::Field<int> pEqn_l_addr = p_rghEqn.lowerAddr();

This does not work:

pEqn_m1.H:50:49: Fehler: »Foam::fvScalarMatrix« hat kein Element namens »upperAddr«
pEqn_m1.H:51:49: Fehler: »Foam::fvScalarMatrix« hat kein Element namens »lowerAddr«

This is weird, since this is what is described in the documentation. I would have understood that the field type is not ok, but not existing?

To get the A matrix, I applied the following:

        double A_pEgn [pEqn_d.size()][pEqn_d.size()];

        forAll(pEqn_d,cellI)
        {
            A_pEgn[cellI][cellI] = pEqn_d[cellI];
        }
        
        forAll(pEqn_u,faceI)
        {
            A_pEgn[pEqn_l_addr[faceI]][pEqn_u_addr[faceI]] = pEqn_u[faceI];
            A_pEgn[pEqn_u_addr[faceI]][pEqn_l_addr[faceI]] = pEqn_l[faceI];
        }

Here at least no complaint. Is there anybody interested in this topic? I only find very old threads with very little code.

Regards,

Daniel