[phuchuynh]

Hi,
I have used power law profile for velocity inlet,
-----------------------------------------------------------

INLET
{
type powerLawTurbulentVelocity;
n (1 0 0);
y (0 1 0);
maxValue 0.032;
alpha 0.224;
zref 0.2;
value uniform (0 0 0);
}

--------------------------------------------------------------
and my decomposeParDict
-----
numberOfSubdomains 2;
method scotch;
-----
I compiled "decomposePar", it Ok, generated two processor*,
when I run parallel " mpirun -np 2 pimpleFoam -parallel ", and got stuck , so as picture below.
Please, help me to solve issue. Thanks

Phuc

[nimasam]

1- look into U file in 0 folder in decomposePar folder and investigate whether zref is there or not

2-decompose the domain in the direction which the whole inlet patch is assigned to only one processor

[phuchuynh]

Originally Posted by nimasam

1- look into U file in 0 folder in decomposePar folder and investigate whether zref is there or not

2-decompose the domain in the direction which the whole inlet patch is assigned to only one processor

Thanks Nimasam,

I have seen the U file in 0 folder in the both processor0, processor1 and I have repaired them by added the zref. So, I compiled "decomposerPar" and run parallel is OK. However, when I use "reconstructPar" , it also have some trouble, and
Could you tell me about it ? How could repair them?

Thanks

Phuc

[nimasam]

well, it seems its not printing zref in U files , maybe you want to add it manually, but it would be time consuming

[phuchuynh]

Originally Posted by nimasam

well, it seems its not printing zref in U files , maybe you want to add it manually, but it would be time consuming

Hi,

Please, I have check and received the success when I run decomposerPar. But I still got stuck in reconstructPar. Could you share me some a ways to adjust it ?

thanks

[wyldckat]

Quick answer @phuchuynh: There are at least 2 ways to solve this problem:

1. Correct the source code for "powerLawTurbulentVelocity". Since you haven't provided the source code for it, it's not possible to tell you right away where it needs to be fixed. Nonetheless, have a look here for ideas: https://github.com/OpenFOAM/OpenFOAM...ryLayer.C#L158
2. Use changeDictionary for repairing the files.

For the second solution, create the file "system/changeDictionaryDict" with the following content:

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      changeDictionaryDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dictionaryReplacement
{
    U
    {
        boundaryField
        {
            INLET
            {
                zref 0.2;
            }
        }
    }
}

Then run the utility like this:

foamJob -p changeDictionary -latestTime

Now you should be able to use reconstructPar like this:

reconstructPar -latestTime

[phuchuynh]

Thank you very much,

This is a my profile, which was changed from parabolicVelocity in OpenFOAM. Please, find a attached file.

// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

//--------------------------------------------------------------------
// Mean velocity and kinetic energy in all faces on inlet patch
//--------------------------------------------------------------------

void powerLawTurbulentVelocityFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    // Get range of inlet patch and face coordinates, c
    boundBox bb(patch().patch().localPoints(), true);
    
    const vectorField& c = patch().Cf();

    // Global y coordinate at face centers of inlet patch
    scalarField coord = (c & y_);

    // Power law profile at inlet with 0 value for v and w components
    vectorField Uinf = (n_*maxValue_*pow (coord, alpha_));
     scalarField Uinfx = Uinf.component(0); //get u component of vector field U

    // Roughness length z0 and turbulent kinetic energy k
    scalar z0 = zref_/exp(1/alpha_);        
    scalarField Iu = 0.1*pow(coord/z0,(-alpha_-0.05));

    scalarField k = 1.5*(pow(Uinfx,2)*pow(Iu,2));
    
//----------------------------------------------------------------------          

    // Calculate local 1-D coordinate for the powerLawTurbulent profile


    vectorField::operator=(Uinf);
}

thanks again.

Phuc

[wyldckat]

Quick answer: OK, if we look at the following constructor:

powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchVectorField(p, iF),
    maxValue_(readScalar(dict.lookup("maxValue"))),
    alpha_(readScalar(dict.lookup("alpha"))),//added line
    zref_(readScalar(dict.lookup("zref"))), //added line

    n_(dict.lookup("n")),
    y_(dict.lookup("y"))
{
    if (mag(n_) < SMALL || mag(y_) < SMALL)
    {
        FatalErrorIn("powerLawTurbulentVelocityFvPatchVectorField(dict)")
            << "n or y given with zero size not correct"
            << abort(FatalError);
    }

    n_ /= mag(n_);
    y_ /= mag(y_);

    evaluate();
}

the following parameters are read from the field file:

    maxValue_(readScalar(dict.lookup("maxValue"))),
    alpha_(readScalar(dict.lookup("alpha"))),//added line
    zref_(readScalar(dict.lookup("zref"))), //added line

    n_(dict.lookup("n")),
    y_(dict.lookup("y"))

If we then look at the code at the end:

// Write
void powerLawTurbulentVelocityFvPatchVectorField::write(Ostream& os) const
{
    fvPatchVectorField::write(os);
    os.writeKeyword("maxValue")
        << maxValue_ << token::END_STATEMENT << nl;

    os.writeKeyword("alpha")
        << alpha_ << token::END_STATEMENT << nl;

    os.writeKeyword("n")
        << n_ << token::END_STATEMENT << nl;
    os.writeKeyword("y")
        << y_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}

The entry "zref_" is missing. It should be changed to something like this:

// Write
void powerLawTurbulentVelocityFvPatchVectorField::write(Ostream& os) const
{
    fvPatchVectorField::write(os);
    os.writeKeyword("maxValue")
        << maxValue_ << token::END_STATEMENT << nl;

    os.writeKeyword("alpha")
        << alpha_ << token::END_STATEMENT << nl;

    os.writeKeyword("zref")
        << zref_ << token::END_STATEMENT << nl;

    os.writeKeyword("n")
        << n_ << token::END_STATEMENT << nl;
    os.writeKeyword("y")
        << y_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}

edit: Don't forget to build after changing and saving the file!!! Namely to build the library with this command:

wmake libso

[hiuluom]

Hi Bruno,

I check the power law velocity with pitzdaily pisoFoam was good, but when I have been running chtmultiregionSimpleFoam in tutorial the boundary met error.

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3   in "/lib/x86_64-linux-gnu/libm.so.6"
#4  pow in "/lib/x86_64-linux-gnu/libm.so.6"
#5  Foam::pow(Foam::Field<double>&, Foam::UList<double> const&, double const&) at ??:?
#6  Foam::pow(Foam::UList<double> const&, double const&) at ??:?
#7  Foam::powerLawTurbulentVelocityFvPatchVectorField::updateCoeffs() at ??:?
#8  Foam::powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField(Foam::fvPatch const&, Foam::DimensionedField<Foam::Vector<double>, Foam::volMesh> const&, Foam::dictionary const&) at ??:?
#9  Foam::fvPatchField<Foam::Vector<double> >::adddictionaryConstructorToTable<Foam::powerLawTurbulentVelocityFvPatchVectorField>::New(Foam::fvPatch const&, Foam::DimensionedField<Foam::Vector<double>, Foam::volMesh> const&, Foam::dictionary const&) at ??:?
#10  Foam::fvPatchField<Foam::Vector<double> >::New(Foam::fvPatch const&, Foam::DimensionedField<Foam::Vector<double>, Foam::volMesh> const&, Foam::dictionary const&) at ??:?
#11  Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::GeometricBoundaryField(Foam::fvBoundaryMesh const&, Foam::DimensionedField<Foam::Vector<double>, Foam::volMesh> const&, Foam::dictionary const&) at ??:?
#12  Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>::readField(Foam::dictionary const&) at ??:?
#13  Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>::readField(Foam::Istream&) at ??:?
#14  Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>::GeometricField(Foam::IOobject const&, Foam::fvMesh const&) at ??:?
#15
 at ??:?
#16  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#17
 at ??:?
Floating point exception (core dumped)

I think it cannot read the parameter in U field. Here is the code in powerLaw

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

\*---------------------------------------------------------------------------*/

#include "powerLawTurbulentVelocityFvPatchVectorField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF
)
:
    fixedValueFvPatchVectorField(p, iF),
    a_(0),
    alpha_(0),
    zref_(0),
    N_(0),
    n_(pTraits<vector>::zero),
    y_(pTraits<vector>::zero)
{}


powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const powerLawTurbulentVelocityFvPatchVectorField& ptf,
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
    a_(ptf.a_),
    alpha_(ptf.alpha_),
    zref_(ptf.zref_),
    N_(ptf.N_),
    n_(ptf.n_),
    y_(ptf.y_)
{}


powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchVectorField(p, iF),
    //maxValue_(readScalar(dict.lookup("maxValue"))),
    a_(readScalar(dict.lookup("a"))),
    alpha_(readScalar(dict.lookup("alpha"))),
    zref_(readScalar(dict.lookup("zref"))),
    N_(readScalar(dict.lookup("N"))),
    n_(dict.lookup("n")),
    y_(dict.lookup("y"))
{
    if (mag(n_) < SMALL || mag(y_) < SMALL)
    {
        FatalErrorIn("powerLawTurbulentVelocityFvPatchVectorField(dict)")
            << "n or y given with zero size not correct"
            << abort(FatalError);
    }

    n_ /= mag(n_);
    y_ /= mag(y_);

    evaluate();
}


powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const powerLawTurbulentVelocityFvPatchVectorField& fcvpvf,
    const DimensionedField<vector, volMesh>& iF
)
:
    fixedValueFvPatchVectorField(fcvpvf, iF),
    //maxValue_(fcvpvf.maxValue_),
    a_(fcvpvf.a_),
    alpha_(fcvpvf.alpha_),
    zref_(fcvpvf.zref_),
    N_(fcvpvf.N_),
    n_(fcvpvf.n_),
    y_(fcvpvf.y_)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void powerLawTurbulentVelocityFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }
    //calculate the global y coordinates of the inlet patch face centers
    boundBox bb(patch().patch().localPoints(), true);
    const vectorField& c = patch().Cf();
    scalarField coord = (c & y_);
    
    vectorField Uinf = (n_*a_*pow(coord,alpha_));
    
    //Turbulent kinetic energy
    
    scalar z0 = zref_/exp(1/alpha_);
    
    scalarField Iu = 1/log(coord/z0);
    
    scalarField Uinfx = Uinf.component(0);
    
    scalarField k = 1.5*(pow(Uinfx,2)*pow(Iu,2));
    
    //calculating the tangential fluctuations
    if (N_ >= c.size())
    {
        Info << "Error: N should be less than or equal to " << c.size() <<
        endl;
    }
    scalarField Af = patch().magSf();
    
    //inlet patch area
    scalar Ap = sum(Af);
    vector bbMax = bb.max();
    scalar lx = bbMax[1];
    scalar ly = bbMax[2];
    scalar L = (2*lx*ly)/(lx+ly);
    
    //vortex size
    scalar sigma = (0.1*L)/2;
    scalar sigma2 = 2*sqr(sigma);    
    
    const scalar t = this->db().time().timeOutputValue();
    if (t >= t)
    {
        vectorField ux = c*0;
        forAll(patch(),facei)
        {
            vector x1 = c[facei];
            scalar x1x = x1[1];
            scalar x1y = (x1 & y_);
            srand (time(NULL) ); vector ux_sum = vector(0, 0, 0);
            for (int iter = 1; iter <= N_; iter ++)
            {
                int RandIndex = rand() % c.size();
                vector xi = c[RandIndex];
                if (xi == x1)
                {
                    scalar Lxi = sqrt(patch().magSf()[iter]);
                    scalar addxi = 0.1*Lxi;
                    xi[1] = xi[1]+addxi;
                    xi[2] = xi[2]+addxi;
                }
                scalar xix = xi[1];
                scalar xiy = (xi & y_);
                scalar dist = sqr(x1x-xix)+sqr(x1y-xiy);
                vector xx = xi-x1;
                scalar ks = k[RandIndex];
                scalar gamma = 4*sqrt((Foam::constant::mathematical::pi*Ap*ks)/(3*N_*0.117783035656383));
                scalar delta = sqrt(patch().magSf()[RandIndex]);
                if (sigma <= delta)
                {
                sigma = delta;
                sigma2 = 2*sqr(sigma);
                }
                vector z1 = vector(1,0,0);
                vector ux1 = (gamma*((xx ^ z1)/dist)*(1-exp(-
                dist/sigma2))*exp(-dist/sigma2));
                ux_sum += ux1;
            }
        ux[facei] = (1/(2*3.141592653589793))*ux_sum;
        }
        vectorField::operator=(Uinf+ux);
    }
}
// Write
void powerLawTurbulentVelocityFvPatchVectorField::write(Ostream& os) const
{
    fvPatchVectorField::write(os);
    // os.writeKeyword("maxValue")
        // << maxValue_ << token::END_STATEMENT << nl;
    os.writeKeyword("a")
        << a_ << token::END_STATEMENT << nl;
    os.writeKeyword("alpha")
        << alpha_ << token::END_STATEMENT << nl;
    os.writeKeyword("zref")
        << zref_ << token::END_STATEMENT << nl;
    os.writeKeyword("N")
        << N_ << token::END_STATEMENT << nl;
    os.writeKeyword("n")
        << n_ << token::END_STATEMENT << nl;
    os.writeKeyword("y")
        << y_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

makePatchTypeField(fvPatchVectorField, powerLawTurbulentVelocityFvPatchVectorField);

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //

How to fixed the error in chtMultiReionSimpleFoam?

Best regards,
Thanh

[hiuluom]

Greeting all,

I modified the code and it can run for multi region or a region.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

\*---------------------------------------------------------------------------*/

#include "powerLawTurbulentVelocityFvPatchVectorField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF
)
:
    fixedValueFvPatchVectorField(p, iF),
    //maxValue_(0),
	Uref_(0),
	alpha_(0),
	zref_(0),
    n_(1, 0, 0),
    y_(0, 1, 0)
{}


powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const powerLawTurbulentVelocityFvPatchVectorField& ptf,
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
    //maxValue_(ptf.maxValue_),
	Uref_(ptf.Uref_),
	alpha_(ptf.alpha_),
	zref_(ptf.zref_),
    n_(ptf.n_),
    y_(ptf.y_)
{}


powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchVectorField(p, iF),
    //maxValue_(readScalar(dict.lookup("maxValue"))),
	Uref_(readScalar(dict.lookup("Uref"))),
	alpha_(readScalar(dict.lookup("alpha"))),
	zref_(readScalar(dict.lookup("zref"))),
    n_(dict.lookup("n")),
    y_(dict.lookup("y"))
{
    if (mag(n_) < SMALL || mag(y_) < SMALL)
    {
        FatalErrorIn("powerLawTurbulentVelocityFvPatchVectorField(dict)")
            << "n or y given with zero size not correct"
            << abort(FatalError);
    }

    n_ /= mag(n_);
    y_ /= mag(y_);

    evaluate();
}


powerLawTurbulentVelocityFvPatchVectorField::powerLawTurbulentVelocityFvPatchVectorField
(
    const powerLawTurbulentVelocityFvPatchVectorField& fcvpvf,
    const DimensionedField<vector, volMesh>& iF
)
:
    fixedValueFvPatchVectorField(fcvpvf, iF),
    //maxValue_(fcvpvf.maxValue_),
	Uref_(fcvpvf.Uref_),
	alpha_(fcvpvf.alpha_),
	zref_(fcvpvf.zref_),
    n_(fcvpvf.n_),
    y_(fcvpvf.y_)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void powerLawTurbulentVelocityFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    // Get range and orientation
    boundBox bb(patch().patch().localPoints(), true);

    //vector ctr = 0.5*(bb.max() + bb.min());
	vector ctr = (bb.min());

    const vectorField& c = patch().Cf();

    // Calculate local 1-D coordinate for the powerLawTurbulent profile
    //scalarField coord = 2*((c - ctr) & y_)/((bb.max() - bb.min()) & y_);
	scalarField coord =((c - ctr) & y_)/((bb.max() - bb.min()) & y_);
	//const scalarField coord = (c & y_);
	
    //vectorField::operator=(n_*Uref_*pow(coord,alpha_));
	vectorField::operator=(n_*Uref_*pow(coord/zref_,alpha_));
	//vectorField Uinf = (n_*Uref_*pow(coord,alpha_));
}


// Write
void powerLawTurbulentVelocityFvPatchVectorField::write(Ostream& os) const
{
    fvPatchVectorField::write(os);
    // os.writeKeyword("maxValue")
        // << maxValue_ << token::END_STATEMENT << nl;
	os.writeKeyword("Uref")
        << Uref_ << token::END_STATEMENT << nl;
	os.writeKeyword("alpha")
        << alpha_ << token::END_STATEMENT << nl;
	os.writeKeyword("zref")
        << zref_ << token::END_STATEMENT << nl;
    os.writeKeyword("n")
        << n_ << token::END_STATEMENT << nl;
    os.writeKeyword("y")
        << y_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

makePatchTypeField(fvPatchVectorField, powerLawTurbulentVelocityFvPatchVectorField);

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //

the result was shown in the figure enclose.

Best,
Thanh

[wyldckat]

Quick answer: You modified a lot the calculations, making it a lot simpler to look at the source code. That certainly makes it easier to diagnose problems in the future.

As for the previous problem, one of the "pow" operations resulted in an invalid operation. I don't which one it was, but you could isolate using old school debugging, i.e. adding messages before and after a calculation is done. For example:

    Info << "Got here 02" << endl;
    Info << "coord: " << coord << endl;
    Info << "alpha_: " << alpha_ << endl;

    vectorField Uinf = (n_*a_*pow(coord,alpha_));

    Info << "Got here 03" << endl;

[rudolf.hellmuth]

Hi, I am having the same problem. The solver with the new BC runs fine in parallel, but it is not reconstructing the solved case because of a division by zero (Foam::sigFpe::sigHandler(int)). My new BC is calculated inlet velocity profile for circular pipes, which takes the flow rate as input. It is a basically the flowRateInletVelocityFvPatchVectorField with a velocity profile (like the parabolic profile).

My patch area is defined as:

// Circular patch area
    scalar A = gSum(patch().magSf());

My debugging with reconstructPar is showing that the area is zero at some point:

...

Reconstructing fields for mesh region0

Time 0.1

Reconstructing FV fields

    Reconstructing volScalarFields

        p

    Reconstructing volVectorFields

        U

        U
***patchArea = 0.787201
***n = 6
***patchRadius = 0.500574
***patchArea = 0
***n = 6
#0 Foam::error:printStack(Foam::Ostream&) in "/opt/OpenFOAM/OpenFOAM-v1606+/platforms/linux64ccDPInt32Opt/lib/libOpenFOAM.so"
#1 Foam::sigFpe::sigHandler(int) in "/opt/OpenFOAM/OpenFOAM-v1606+/platforms/linux64ccDPInt32Opt/lib/libOpenFOAM.so"

...

(1) I don't know why it calculating the BC during parallel reconstruction.
(2) Why is it calculating it twice?
(3) Why it is calculating it as A = 0?

Thanks
Rudolf

[rudolf.hellmuth]

I've fixed the issue, but I haven't understood exactly what was going on. Please tell me, if you know.

The fix was to change (original):

Foam::generalTubeVelocityProfileFvPatchVectorField::
generalTubeVelocityProfileFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchVectorField(p, iF),
    n_(dict.lookupOrDefault<scalar>("n", 2.0)),
    Q_(Function1<scalar>::New("Q", dict))
{
        evaluate();
}

To (working version):

Foam::generalTubeVelocityProfileFvPatchVectorField::
generalTubeVelocityProfileFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchVectorField(p, iF),
    n_(dict.lookupOrDefault<scalar>("n", 2.0)),
    Q_(Function1<scalar>::New("Q", dict))
{
    if (dict.found("value"))
    {
        fvPatchField<vector>::operator=
        (
            vectorField("value", dict, p.size())
        );
        Info << "******IF******" << endl;
    }
    else
    {
        evaluate(Pstream::blocking);
        Info << "******ELSE******" << endl;
    }
}

* Notice the debugging tags.

I haven't written any "value" in the BC dictionary. Any application that reads the velocity field first would pop the "******ELSE******" tag. Then, all subsequent applications, would pop the "******IF******" tag. Thus, evaluate(Pstream::blocking); is called once, by any application, and it saves the value nonuniform <vector list> on ./0/U.

When it comes to reconstructPar, now it is calling updateCoeffs just once:

Reconstructing fields for mesh region0

Time = 0

Reconstructing FV fields

    Reconstructing volScalarFields

        p

    Reconstructing volVectorFields

        U
******IF******
******IF******
******IF******
******IF******
******IF******
******IF******
******IF******
******IF******
****A = 0.787201
****n = 6
****R = 0.500574

    Reconstructing surfaceScalarFields

        phi

*Notice that the "******IF******" tag is being called once for each of the 8 processors.


Again, it is reconstructing in parallel now, but I don't know why. I have just copied that fix from the original code. I had thought that that was unnecessary.

Best regards,
Rudolf