CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

overwrite matrix / vector value

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   June 24, 2014, 10:56
Default overwrite matrix / vector value
  #1
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 141
Rep Power: 7
danny123 is on a distinguished road
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:

Code:
        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):

HTML Code:
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 is offline   Reply With Quote

Old   June 30, 2014, 06:01
Default Update
  #2
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 141
Rep Power: 7
danny123 is on a distinguished road
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:

Code:
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:

HTML Code:
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:
Code:
    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 is offline   Reply With Quote

Old   July 1, 2014, 05:32
Default A and B in FvScalarMatrix
  #3
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 141
Rep Power: 7
danny123 is on a distinguished road
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:

Code:
        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:

Code:
        Foam::Field<int> pEqn_u_addr = p_rghEqn.upperAddr();
        Foam::Field<int> pEqn_l_addr = p_rghEqn.lowerAddr();
This does not work:

HTML Code:
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:

Code:
        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
danny123 is offline   Reply With Quote

Reply

Tags
matrix value overwrite

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
unknown asymmetric matrix solver ebaldwin OpenFOAM Running, Solving & CFD 2 February 17, 2014 14:07
more equation in block matrix system yhaomin2007 OpenFOAM 1 September 6, 2012 08:33
Force can not converge colopolo CFX 13 October 4, 2011 22:03
OpenFOAM version 1.6 details lakeat OpenFOAM Running, Solving & CFD 42 August 26, 2009 21:47
Elemtary matrix to CSR global matrix xueying Main CFD Forum 2 September 24, 2002 09:44


All times are GMT -4. The time now is 11:45.