Non classé

OpenFOAM – fixedTabulatedVectorValue

 

This boundary condition allows to impose vector values stored in a table on a side of the fluid domain. Vector values must be associated with a point, defined by its coordinates.

x y z Ux Uy Uz
-12 -42 0 -0.0124975 -0.456742 0.00380209
-12 -42 0.274078 -0.0248316 -0.727996 0.00482216
-12 -41.6064 0.274078 -0.011896 -0.45603 0.00395441

The program uses a closest neighbours algorithm to associate an input value calculated by the inverse distance weighting method to each face center of the patch on which the condition is imposed

The developments use the nanoflann library accessible by clicking on the following link. This library does not require any previous compilation, however the links to the headers must be correctly defined in order to compile the sources of fixedTabularVectorValue

In the following, we provide the sources of the .h and .cpp files as well as the two option files and files needed for compilation. The paths are likely to change depending on the work environment.

 

  • options
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(FOAM_DEV)/libraries/nanoflann/include \
-I$(FOAM_DEV)/libraries/utils/include
EXE_LIBS =<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • files
fvPatchFields = ./
derivedFvPatchFields = $(fvPatchFields)/derived

$(derivedFvPatchFields)/fixedTabulatedVectorValue/fixedTabulatedVectorValueFvPatchScalarField.C

LIB = $(FOAM_USER_LIBBIN)/libArepCustomPatch<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • fixedTabulatedVectorValueFvPatchScalarField.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
 * Copyright 2018 arep <alexis.sauvageon@arep.fr>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above
 *   copyright notice, this list of conditions and the following disclaimer
 *   in the documentation and/or other materials provided with the
 *   distribution.
 * * Neither the name of the  nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Class
    fixedTabulatedVectorValueFvPatchVectorField

Description
    Compute boundaries conditions from a field stored in table (format .raw)
	[X Y Z VX VY VZ]
    the two first lines of the table are not interpreted.
    The interpollated value is computed thanks to a PID method.

    \heading Patch usage

    \table
        Property    | Description             		| Required   	| Default value
        fileDir     | path to the input file  		| yes        	| 
        nbPoint     | number of point for PID 	    | no         	| 4
        wFactor     | p factor for PID      		| no        	| 2.0

    \endtable

    Example of the boundary condition specification:
    \verbatim
    myPatch
    {
        type            fixedTabulatedVectorValue;
        fileDir         "constant/inputTableData/p_boxValue.raw";
        nbPoint              8;
        wFactor            2.0;

    }
    \endverbatim


SourceFiles
    fixedTabulatedVectorValueFvPatchVectorField.C

Author
    Alexis Sauvageon.  All rights reserved

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

#ifndef fixedTabulatedVectorValueFvPatchVectorField_H
#define fixedTabulatedVectorValueFvPatchVectorField_H

#include "fvPatchFields.H"
#include "fixedValueFvPatchFields.H"

#include <fstream>
#include <sstream>
#include <iterator>



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

namespace Foam
{

/*---------------------------------------------------------------------------*\
              Class fixedTabulatedVectorValueFvPatchField Declaration
\*---------------------------------------------------------------------------*/

class fixedTabulatedVectorValueFvPatchVectorField
:
    public fixedValueFvPatchVectorField
{
    // Private data

        //- Peak velocity magnitude
        std::string fileDir_;
	unsigned int nbPoint_;
	scalar wFactor_;

public:

    //- Runtime type information
    TypeName("fixedTabulatedVectorValue");


    // Constructors

        //- Construct from patch and internal field
        fixedTabulatedVectorValueFvPatchVectorField
        (
            const fvPatch&,
            const DimensionedField<vector, volMesh>&
        );

        //- Construct from patch, internal field and dictionary
        fixedTabulatedVectorValueFvPatchVectorField
        (
            const fvPatch&,
            const DimensionedField<vector, volMesh>&,
            const dictionary&
        );

        //- Construct by mapping given fixedTabulatedVectorValueFvPatchVectorField
        //  onto a new patch
        fixedTabulatedVectorValueFvPatchVectorField
        (
            const fixedTabulatedVectorValueFvPatchVectorField&,
            const fvPatch&,
            const DimensionedField<vector, volMesh>&,
            const fvPatchFieldMapper&
        );

        //- Construct and return a clone
        virtual tmp<fvPatchVectorField> clone() const
        {
            return tmp<fvPatchVectorField>
            (
                new fixedTabulatedVectorValueFvPatchVectorField(*this)
            );
        }

        //- Construct as copy setting internal field reference
        fixedTabulatedVectorValueFvPatchVectorField
        (
            const fixedTabulatedVectorValueFvPatchVectorField&,
            const DimensionedField<vector, volMesh>&
        );

        //- Construct and return a clone setting internal field reference
        virtual tmp<fvPatchVectorField> clone
        (
            const DimensionedField<vector, volMesh>& iF
        ) const
        {
            return tmp<fvPatchVectorField>
            (
                new fixedTabulatedVectorValueFvPatchVectorField(*this, iF)
            );
        }


    // Member functions

        //- Return fileDir_
        std::string& fileDir()
        {
            return fileDir_;
        }
        //- Return nbPoint_
        unsigned int& nbPoint()
        {
            return nbPoint_;
        }
        //- Return wFactor_
        scalar& wFactor()
        {
            return wFactor_;
        }


        //- Update coefficients
        virtual void updateCoeffs();

        //- Write
        virtual void write(Ostream&) const;
};


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

} // End namespace Foam

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

#endif

// ************************************************************************* //
<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span><span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • fixedTabulatedVectorValueFvPatchScalarField.C
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
 * Copyright 2018 arep <alexis.sauvageon@arep.fr>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above
 *   copyright notice, this list of conditions and the following disclaimer
 *   in the documentation and/or other materials provided with the
 *   distribution.
 * * Neither the name of the  nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

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

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


#include "nanoflann.hpp"
#include "kdtree.h"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


namespace Foam
{

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

fixedTabulatedVectorValueFvPatchVectorField::fixedTabulatedVectorValueFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF
)
:
    fixedValueFvPatchVectorField(p, iF),
    fileDir_("."),
    nbPoint_(4),
    wFactor_(2)
{}


fixedTabulatedVectorValueFvPatchVectorField::fixedTabulatedVectorValueFvPatchVectorField
(
    const fixedTabulatedVectorValueFvPatchVectorField& ptf,
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedValueFvPatchVectorField(ptf, p, iF, mapper),
    fileDir_(ptf.fileDir_),
    nbPoint_(ptf.nbPoint_),
    wFactor_(ptf.wFactor_)
{}


fixedTabulatedVectorValueFvPatchVectorField::fixedTabulatedVectorValueFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchVectorField(p, iF),
    fileDir_(string(dict.lookup("fileDir"))),
    nbPoint_(dict.lookupOrDefault<scalar>("nbPoint",4)),
    wFactor_(dict.lookupOrDefault<scalar>("wFactor",2))
{

// check the existance of the file

    std::ifstream inFile(fileDir_,std::ios::in);  
 
    if(inFile) 
    {       
        inFile.close(); 
    }
    else 
    {
        FatalErrorIn("fixedTabulatedVectorValueFvPatchVectorField(dict)")
            << "path to the raw file is not correct"
            << abort(FatalError);
    }

    if( (nbPoint_ < 1) || (wFactor_ <= 0) ) 
    {
        FatalErrorIn("fixedTabulatedVectorValueFvPatchVectorField(dict)")
            << "the number of interpolation point is less than 1 or the weighting factor is not positive."
            << abort(FatalError);
    }

    evaluate();
}


fixedTabulatedVectorValueFvPatchVectorField::fixedTabulatedVectorValueFvPatchVectorField
(
    const fixedTabulatedVectorValueFvPatchVectorField& fcvpvf,
    const DimensionedField<vector, volMesh>& iF
)
:
    fixedValueFvPatchVectorField(fcvpvf, iF),
    fileDir_(fcvpvf.fileDir_),
    nbPoint_(fcvpvf.nbPoint_),
    wFactor_(fcvpvf.wFactor_)
{}


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

void fixedTabulatedVectorValueFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }
// >> open the input data file
	std::ifstream inFile(fileDir_,std::ios::in);
        // count the lines in the file
	inFile.unsetf(std::ios_base::skipws);
    	unsigned nbLines = std::count(
        	std::istream_iterator<char>(inFile),
        	std::istream_iterator<char>(), 
        	'\n');
	nbLines=nbLines-2;
	//...skip two first lines (commentary)
	inFile.close();
	inFile.open(fileDir_);
	std::string line;  
        std::getline(inFile, line);  
        std::getline(inFile, line);  

	//Initialize a plain continous array for the data
        scalar vectorData[nbLines][6];
	unsigned int index_(0);

	PointCloud<scalar> cloud;
	//...for all line until the end, create a VectorField with the computed point
        while(std::getline(inFile, line))  // tant que l'on peut mettre la ligne dans "contenu"
        {
		std::istringstream buffer(line);
		std::string word_("");

	        buffer >> word_;
		scalar x_(atof(word_.c_str()));
	        buffer >> word_;
		scalar y_(atof(word_.c_str()));
	        buffer >> word_;
		scalar z_(atof(word_.c_str()));
	        buffer >> word_;
		scalar valx_(atof(word_.c_str()));
	        buffer >> word_;
		scalar valy_(atof(word_.c_str()));
	        buffer >> word_;
		scalar valz_(atof(word_.c_str()));

		updatePointCloud(cloud, x_, y_, z_); 

		vectorData[index_][0] = x_;
                vectorData[index_][1] = y_;
                vectorData[index_][2] = z_;
                vectorData[index_][3] = valx_;
		vectorData[index_][4] = valy_;
		vectorData[index_][5] = valz_;
		index_++;
	
        }

	
	// construct a kd-tree index:
	typedef nanoflann::KDTreeSingleIndexAdaptor<
		nanoflann::L2_Simple_Adaptor<scalar, PointCloud<scalar> > ,
		PointCloud<scalar>,
		3
		> my_kd_tree_t;

	my_kd_tree_t   index(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(10) );
	index.buildIndex();

#if 0
	// Test resize of dataset and rebuild of index:
	cloud.pts.resize(cloud.pts.size()*0.5);
	index.buildIndex();
#endif


// >> get the face of the current patch & initialize a VectorField of the right size
	const vectorField& faceCenter = patch().Cf();
	vectorField values(faceCenter.size(), vector(0,0,0));

// >> compute the scalar value for all face centers
        forAll(faceCenter, f)
        {
// >>>>>>> get face center point
                const point faceCenterPoint = faceCenter[f];
// >>>>>>> find closest neighbourg
		// define the requested point
		const scalar query_pt[3] = { faceCenterPoint.x(), faceCenterPoint.y(), faceCenterPoint.z()};
		// find the 8 closest points
		size_t num_results = nbPoint_;
		std::vector<size_t>   ret_index(num_results);
		std::vector<scalar> out_dist_sqr(num_results);

		num_results = index.knnSearch(&query_pt[0], num_results, &ret_index[0], &out_dist_sqr[0]);
		
		// In case of less points in the tree than requested:
		ret_index.resize(num_results);
		out_dist_sqr.resize(num_results);

		// weight factors for PID
		std::vector<scalar> out_dist_weight(num_results);
		for (size_t i = 0; i < num_results; i++)
		{
			out_dist_weight[i]=1/std::pow(out_dist_sqr[i],wFactor_);
		}	
		// Distance inverse weighting methode (PID)
		scalar WxScalar(0);
		scalar WyScalar(0);
		scalar WzScalar(0);
		scalar Wtot(0);
		for (size_t i = 0; i < num_results; i++)
		{
			WxScalar+=out_dist_weight[i]*vectorData[ret_index[i]][3];
			WyScalar+=out_dist_weight[i]*vectorData[ret_index[i]][4];
			WzScalar+=out_dist_weight[i]*vectorData[ret_index[i]][5];
			Wtot+=out_dist_weight[i];
		}
		values[f]=vector(WxScalar/Wtot, WyScalar/Wtot, WzScalar/Wtot); ;


        }



    // Calculate local 1-D coordinate for the parabolic profile
   // VectorField coord = 0; //2*((c - ctr) & y_)/((bb.max() - bb.min()) & y_);

   vectorField::operator=(values);
}


// Write
void fixedTabulatedVectorValueFvPatchVectorField::write(Ostream& os) const
{
    fvPatchVectorField::write(os);
    os.writeKeyword("fileDir")
        << fileDir_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}


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

makePatchTypeField(fvPatchVectorField, fixedTabulatedVectorValueFvPatchVectorField);

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

} // End namespace Foam

// ************************************************************************* //
<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span><span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

 

Non classé

OpenFOAM – fixedTabulatedScalarValue

This boundary condition allows to impose vector values stored in a table on a side of the fluid domain. Scalar values must be associated with a point, defined by its coordinates.

x y z epsilon
-12 -42 0 0.00800494
-12 -42 0.274078 0.00730994
-12 -41.6064 0.274078 0.00792009

The program uses a closest neighbours algorithm to associate an input value calculated by the inverse distance weighting method to each face center of the patch on which the condition is imposed

The developments use the nanoflann library accessible by clicking on the following link. This library does not require any previous compilation, however the links to the headers must be correctly defined in order to compile the sources of fixedTabularVectorValue

In the following, we provide the sources of the .h and .cpp files as well as the two option files and files needed for compilation. The paths are likely to change depending on the work environment.

 

  • options
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(FOAM_DEV)/libraries/nanoflann/include \
-I$(FOAM_DEV)/libraries/utils/include
EXE_LIBS =<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • files
fvPatchFields = ./
derivedFvPatchFields = $(fvPatchFields)/derived

$(derivedFvPatchFields)/fixedTabulatedScalarValue/fixedTabulatedScalarValueFvPatchScalarField.C

LIB = $(FOAM_USER_LIBBIN)/libArepCustomPatch

 

  • fixedTabulatedScalarValueFvPatchScalarField.H
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
 * Copyright 2018 arep <alexis.sauvageon@arep.fr>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above
 *   copyright notice, this list of conditions and the following disclaimer
 *   in the documentation and/or other materials provided with the
 *   distribution.
 * * Neither the name of the  nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Class
    fixedTabulatedScalarValueFvPatchScalarField

Description
    Compute boundaries conditions from a field stored in table (format .raw)
	[X Y Z Scalar]
    the two first lines of the table are not interpreted.
    The interpollated value is computed thanks to a PID method.

    \heading Patch usage

    \table
        Property    | Description             	 	| Required   	| Default value
        fileDir     | path to the input file  		| yes        	| 
        nbPoint     | number of point for PID 	        | no       	| 4
        wFactor     | p factor for PID      		| no        	| 2.0

    \endtable

    Example of the boundary condition specification:
    \verbatim
    myPatch
    {
        type            fixedTabulatedScalarValue;
        fileDir         "constant/inputTableData/p_boxValue.raw";
        nbPoint              8;
        wFactor            2.0;

    }
    \endverbatim


SourceFiles
    fixedTabulatedScalarValueFvPatchScalarField.C

Author
    Alexis Sauvageon.  All rights reserved

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

#ifndef fixedTabulatedScalarValueFvPatchScalarField_H
#define fixedTabulatedScalarValueFvPatchScalarField_H

#include "fvPatchFields.H"
#include "fixedValueFvPatchFields.H"


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

namespace Foam
{

/*---------------------------------------------------------------------------*\
              Class fixedTabulatedScalarValueFvPatchField Declaration
\*---------------------------------------------------------------------------*/

class fixedTabulatedScalarValueFvPatchScalarField
:
    public fixedValueFvPatchScalarField
{
    // Private data

        //- Peak velocity magnitude
        std::string fileDir_;
	unsigned int nbPoint_;
	scalar wFactor_;

public:

    //- Runtime type information
    TypeName("fixedTabulatedScalarValue");


    // Constructors

        //- Construct from patch and internal field
        fixedTabulatedScalarValueFvPatchScalarField
        (
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&
        );

        //- Construct from patch, internal field and dictionary
        fixedTabulatedScalarValueFvPatchScalarField
        (
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&,
            const dictionary&
        );

        //- Construct by mapping given fixedTabulatedScalarValueFvPatchScalarField
        //  onto a new patch
        fixedTabulatedScalarValueFvPatchScalarField
        (
            const fixedTabulatedScalarValueFvPatchScalarField&,
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&,
            const fvPatchFieldMapper&
        );

        //- Construct and return a clone
        virtual tmp<fvPatchScalarField> clone() const
        {
            return tmp<fvPatchScalarField>
            (
                new fixedTabulatedScalarValueFvPatchScalarField(*this)
            );
        }

        //- Construct as copy setting internal field reference
        fixedTabulatedScalarValueFvPatchScalarField
        (
            const fixedTabulatedScalarValueFvPatchScalarField&,
            const DimensionedField<scalar, volMesh>&
        );

        //- Construct and return a clone setting internal field reference
        virtual tmp<fvPatchScalarField> clone
        (
            const DimensionedField<scalar, volMesh>& iF
        ) const
        {
            return tmp<fvPatchScalarField>
            (
                new fixedTabulatedScalarValueFvPatchScalarField(*this, iF)
            );
        }


    // Member functions

        //- Return fileDir_
        std::string& fileDir()
        {
            return fileDir_;
        }
        //- Return nbPoint_
        unsigned int& nbPoint()
        {
            return nbPoint_;
        }
        //- Return wFactor_
        scalar& wFactor()
        {
            return wFactor_;
        }


        //- Update coefficients
        virtual void updateCoeffs();

        //- Write
        virtual void write(Ostream&) const;
};


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

} // End namespace Foam

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

#endif

// ************************************************************************* //
<span id="mce_marker" data-mce-type="bookmark" data-mce-fragment="1">​</span>

 

  • fixedTabulatedScalarValueFvPatchScalarField.C
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
 * Copyright 2018 arep <alexis.sauvageon@arep.fr>
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 
 * * Redistributions of source code must retain the above copyright
 *   notice, this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above
 *   copyright notice, this list of conditions and the following disclaimer
 *   in the documentation and/or other materials provided with the
 *   distribution.
 * * Neither the name of the  nor the names of its
 *   contributors may be used to endorse or promote products derived from
 *   this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

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

#include "fixedTabulatedScalarValueFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "volFields.H"
#include "surfaceFields.H"
#include <fstream>
#include <sstream>
#include <iterator>

#include "nanoflann.hpp"
#include "kdtree.h"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


namespace Foam
{

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

fixedTabulatedScalarValueFvPatchScalarField::fixedTabulatedScalarValueFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedValueFvPatchScalarField(p, iF),
    fileDir_("."),
    nbPoint_(4),
    wFactor_(2)
{}


fixedTabulatedScalarValueFvPatchScalarField::fixedTabulatedScalarValueFvPatchScalarField
(
    const fixedTabulatedScalarValueFvPatchScalarField& ptf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
    fileDir_(ptf.fileDir_),
    nbPoint_(ptf.nbPoint_),
    wFactor_(ptf.wFactor_)
{}


fixedTabulatedScalarValueFvPatchScalarField::fixedTabulatedScalarValueFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchScalarField(p, iF),
    fileDir_(string(dict.lookup("fileDir"))),
    nbPoint_(dict.lookupOrDefault<scalar>("nbPoint",4)),
    wFactor_(dict.lookupOrDefault<scalar>("wFactor",2))
{

// check the existance of the file

    std::ifstream inFile(fileDir_,std::ios::in);  
 
    if(inFile) 
    {       
        inFile.close(); 
    }
    else 
    {
        FatalErrorIn("fixedTabulatedScalarValueFvPatchScalarField(dict)")
            << "path to the raw file is not correct"
            << abort(FatalError);
    }

    if( (nbPoint_ < 1) || (wFactor_ <= 0) ) 
    {
        FatalErrorIn("fixedTabulatedScalarValueFvPatchScalarField(dict)")
            << "the number of interpolation point is less than 1 or the weighting factor is not positive."
            << abort(FatalError);
    }

    evaluate();
}


fixedTabulatedScalarValueFvPatchScalarField::fixedTabulatedScalarValueFvPatchScalarField
(
    const fixedTabulatedScalarValueFvPatchScalarField& fcvpvf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedValueFvPatchScalarField(fcvpvf, iF),
    fileDir_(fcvpvf.fileDir_),
    nbPoint_(fcvpvf.nbPoint_),
    wFactor_(fcvpvf.wFactor_)
{}


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

void fixedTabulatedScalarValueFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }
// >> open the input data file
	std::ifstream inFile(fileDir_,std::ios::in);
        // count the lines in the file
	inFile.unsetf(std::ios_base::skipws);
    	unsigned nbLines = std::count(
        	std::istream_iterator<char>(inFile),
        	std::istream_iterator<char>(), 
        	'\n');
	nbLines=nbLines-2;
	//...skip two first lines (commentary)
	inFile.close();
	inFile.open(fileDir_);
	std::string line;  
        std::getline(inFile, line);  
        std::getline(inFile, line);  

	//Initialize a plain continous array for the data
        scalar vectorData[nbLines][4];
	unsigned int index_(0);

	PointCloud<scalar> cloud;
	//...for all line until the end, create a scalarField with the computed point
        while(std::getline(inFile, line))  // tant que l'on peut mettre la ligne dans "contenu"
        {
		std::istringstream buffer(line);
		std::string word_("");

	        buffer >> word_;
		scalar x_(atof(word_.c_str()));
	        buffer >> word_;
		scalar y_(atof(word_.c_str()));
	        buffer >> word_;
		scalar z_(atof(word_.c_str()));
	        buffer >> word_;
		scalar val_(atof(word_.c_str()));

		updatePointCloud(cloud, x_, y_, z_); 

		vectorData[index_][0] = x_;
                vectorData[index_][1] = y_;
                vectorData[index_][2] = z_;
                vectorData[index_][3] = val_;
		index_++;
	
        }

	
	// construct a kd-tree index:
	typedef nanoflann::KDTreeSingleIndexAdaptor<
		nanoflann::L2_Simple_Adaptor<scalar, PointCloud<scalar> > ,
		PointCloud<scalar>,
		3
		> my_kd_tree_t;

	my_kd_tree_t   index(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(10) );
	index.buildIndex();

#if 0
	// Test resize of dataset and rebuild of index:
	cloud.pts.resize(cloud.pts.size()*0.5);
	index.buildIndex();
#endif


// >> get the face of the current patch & initialize a scalarField of the right size
	const vectorField& faceCenter = patch().Cf();
	scalarField values(faceCenter.size(), scalar(0));

// >> compute the scalar value for all face centers
        forAll(faceCenter, f)
        {
// >>>>>>> get face center point
                const point faceCenterPoint = faceCenter[f];
// >>>>>>> find closest neighbourg
		// define the requested point
		const scalar query_pt[3] = { faceCenterPoint.x(), faceCenterPoint.y(), faceCenterPoint.z()};
		// find the 8 closest points
		size_t num_results = nbPoint_;
		std::vector<size_t>   ret_index(num_results);
		std::vector<scalar> out_dist_sqr(num_results);

		num_results = index.knnSearch(&query_pt[0], num_results, &ret_index[0], &out_dist_sqr[0]);
		
		// In case of less points in the tree than requested:
		ret_index.resize(num_results);
		out_dist_sqr.resize(num_results);

		// weight factors for PID
		std::vector<scalar> out_dist_weight(num_results);
		for (size_t i = 0; i < num_results; i++)
		{
			out_dist_weight[i]=1/std::pow(out_dist_sqr[i],wFactor_);
		}	
		// Distance inverse weighting methode (PID)
		scalar WxScalar(0);
		scalar Wtot(0);
		for (size_t i = 0; i < num_results; i++)
		{
			WxScalar+=out_dist_weight[i]*vectorData[ret_index[i]][3];
			Wtot+=out_dist_weight[i];
		}
		values[f]=WxScalar/Wtot;

		// print the index of the closest point and their distance
		/*cout << faceCenterPoint.x() << "    " << faceCenterPoint.y() << "    " << faceCenterPoint.z()  << std::endl;
		cout << values[f] << std::endl;
		cout << std::endl;

		for (size_t i = 0; i < num_results; i++)
		{
		cout << vectorData[ret_index[i]][0]  << "    " << vectorData[ret_index[i]][1] << "    " << vectorData[ret_index[i]][2] << std::endl;
		cout << "dist : " << out_dist_sqr[i] << " value : " << vectorData[ret_index[i]][4]  << std::endl;
		cout << "weight : " << out_dist_weight[i]   << std::endl;
		}*/

        }



    // Calculate local 1-D coordinate for the parabolic profile
   // scalarField coord = 0; //2*((c - ctr) & y_)/((bb.max() - bb.min()) & y_);

   scalarField::operator=(values);
}


// Write
void fixedTabulatedScalarValueFvPatchScalarField::write(Ostream& os) const
{
    fvPatchScalarField::write(os);
    os.writeKeyword("fileDir")
        << fileDir_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}


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

makePatchTypeField(fvPatchScalarField, fixedTabulatedScalarValueFvPatchScalarField);

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

} // End namespace Foam

// ************************************************************************* //
Non classé

Calculation of the incident solar flux on a cylinder

And a routine for calculating the incident solar flux on an individual, represented by a vertical cylinder (method described in this article).

###################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
#
#  AREP, 16 avenue d'Ivry,75013 Paris, FRANCE
##################################################################
#Calcul du flux recu par un cylindre
#	calcul integral pour le flux direct
#	modelisation anisotrope pour le flux diffus,
#   d'apres la methode de (Perez et. al 1990)
##################################################################
# Last modification : 01/08/2017                                 #
##################################################################
# Copyright (C) 2018 El Mehdi Hamdani et Edouard Walther         #  
# This program 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   #
# of the License, or (at your option) any later version.         #
##################################################################
# Contact : elmehdi[dot]hamdani[at]arep[dot]fr 
#           edouard[dot]walther[at]arep[dot]fr
##################################################################

import math
global pi
pi = 3.141592653589793

# Cette fonction permet de retourner la valeur du flux diffus circumsolaire, a partir des valeurs du flux diffus horizontal, du flux direct normal, de la hauteur solaire et de l'angle d'incidence
def flux_diffus_circumsolaire_perez(diffus_horizontal, direct_normal, hauteur_solaire, angle_incidence):

	# Valeur de la constante solaire (flux solaire moyen qui parvient sur la terre a la limite de l'atmosphere) en W/m²
	I0 = 1368
	# angle solaire zenithal
	ksi_z = 90 - math.degrees(hauteur_solaire)
	# composante decrivant la transparence du ciel
	epsilon = (((diffus_horizontal + direct_normal)/diffus_horizontal) + 5.535 * 0.000001 * math.pow(ksi_z, 3))/(1 + 5.535 * 0.000001 * math.pow(ksi_z, 3))
	# Masse atmospherique (Kasten and Young 1989)
	a = max(0, math.cos(angle_incidence))
	b = max(math.cos(math.radians(85)), math.cos(math.radians(ksi_z)))
	ma = 1/(math.cos(math.radians(ksi_z)) + 0.50572 * ((96.07995 - ksi_z)**-1.6364))
	

	# Delta represente l'eclairement du ciel
	delta = ma * diffus_horizontal / I0
	
	# coefficients donnes par PEREZ et al (1990)
	if epsilon <= 1.065:
		f11 = 0.013
		f12 = 0.764
		f13 = -0.1
		f21 = -0.058
		f22 = 0.127
		f23 = -0.023
	elif epsilon <= 1.23:
		f11 = 0.095
		f12 = 0.920
		f13 = -0.152
		f21 = 0
		f22 = 0.051
		f23 = -0.02
	elif epsilon <= 1.5:
		f11 = 0.464
		f12 = 0.421
		f13 = -0.28
		f21 = 0.064
		f22 = -0.051
		f23 = -0.002
	elif epsilon <= 1.95:
		f11 = 0.759
		f12 = -0.009
		f13 = -0.373
		f21 = 0.201
		f22 = -0.382
		f23 = 0.01
	elif epsilon <= 2.8:
		f11 = 0.976
		f12 = -0.4
		f13 = -0.436
		f21 = 0.271
		f22 = -0.638
		f23 = 0.051
	elif epsilon <= 4.5:
		f11 = 1.176
		f12 = -1.254
		f13 = -0.462
		f21 = 0.295
		f22 = -0.975
		f23 = 0.129
	elif epsilon <= 6.2:
		f11 = 1.106
		f12 = -1.563
		f13 = -0.398
		f21 = 0.301
		f22 = -1.442
		f23 = 0.212
	else:
		f11 = 0.934
		f12 = -1.501
		f13 = -0.271
		f21 = 0.42
		f22 = -2.917
		f23 = 0.249
	
	# coefficient ponderant les densites du flux circumsolaire
	K1 = max(0, f11 + f12 * delta + pi * ksi_z * f13 / 180)
	
	# coefficient ponderant les densites du flux de l'horizon
	K2 = f21 + f22 * delta + pi * ksi_z * f23 / 180
	
	# flux diffus provenant du rayonnement circumsolaire W/m²
	diffus_circumsolaire = diffus_horizontal * K1 * a/b
	
	#retourne la valeur du flux diffus circumsolaire
	return diffus_circumsolaire

def flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, angle_incidence):
# Cette fonction permet de retourner la valeur du flux diffus circumsolaire, a partir des valeurs du flux diffus horizontal, du flux direct normal, de la hauteur solaire et de l'angle d'incidence

	# Valeur de la constante solaire (flux solaire moyen qui parvient sur la terre a la limite de l'atmosphere) en W/m²
	I0 = 1368
	# angle solaire zenithal
	ksi_z = 90 - math.degrees(hauteur_solaire)
	# composante decrivant la transparence du ciel
	epsilon = (((diffus_horizontal + direct_normal)/diffus_horizontal) + 5.535 * 0.000001 * math.pow(ksi_z, 3))/(1 + 5.535 * 0.000001 * math.pow(ksi_z, 3))
	
	# Masse atmospherique (Kasten and Young 1989)
	a = max(0, math.cos(angle_incidence))
	b = max(math.cos(math.radians(85)), math.cos(math.radians(ksi_z)))
	ma = 1/(math.cos(math.radians(ksi_z)) + 0.50572 * ((96.07995 - ksi_z)**-1.6364))
	
	# Delta represente l'eclairement du ciel
	delta = ma * diffus_horizontal / I0
	
	# coefficients donnes par PEREZ et al (1990)
	if epsilon <= 1.065:
		f11 = 0.013
		f12 = 0.764
		f13 = -0.1
		f21 = -0.058
		f22 = 0.127
		f23 = -0.023
		
	elif epsilon <= 1.23:
		f11 = 0.095
		f12 = 0.920
		f13 = -0.152
		f21 = 0
		f22 = 0.051
		f23 = -0.02
		
	elif epsilon <= 1.5:
		f11 = 0.464
		f12 = 0.421
		f13 = -0.28
		f21 = 0.064
		f22 = -0.051
		f23 = -0.002
		
	elif epsilon <= 1.95:
		f11 = 0.759
		f12 = -0.009
		f13 = -0.373
		f21 = 0.201
		f22 = -0.382
		f23 = 0.01
		
	elif epsilon <= 2.8:
		f11 = 0.976
		f12 = -0.4
		f13 = -0.436
		f21 = 0.271
		f22 = -0.638
		f23 = 0.051
		
	elif epsilon <= 4.5:
		f11 = 1.176
		f12 = -1.254
		f13 = -0.462
		f21 = 0.295
		f22 = -0.975
		f23 = 0.129
		
	elif epsilon <= 6.2:
		f11 = 1.106
		f12 = -1.563
		f13 = -0.398
		f21 = 0.301
		f22 = -1.442
		f23 = 0.212
		
	else:
		f11 = 0.934
		f12 = -1.501
		f13 = -0.271
		f21 = 0.42
		f22 = -2.917
		f23 = 0.249
		
	# coefficient ponderant les densites du flux circumsolaire
	K1 = max(0, f11 + f12 * delta + pi * ksi_z * f13 / 180)
	
	# coefficient ponderant les densites du flux de l'horizon
	K2 = f21 + f22 * delta + pi * ksi_z * f23 / 180
	
	# flux diffus provenant de l'horizon et de la voute celeste W/m²
	diffus_horizon_voute = diffus_horizontal * ((1 - K1) * (0.5 * (1 + math.cos(pi/2))) + K2 * math.sin(pi/2))
	
	# retourne la valeur du flux diffus de l'horizon et de la voute celeste
	return diffus_horizon_voute
	
# Cette fonction permet de retourner la valeur du flux solaire absorbe par un individu, 
# a partir des valeurs des flux solaires diffus et direct horizontaux, de l'albedo du sol, du jour de l'annee, de l'heure, du numero du fuseau, de la latitude et de la longitude du site
def flux_cylindre_anisotrope (diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est): 

	# source : Wikipedia. L'equation du temps permet de calculer le temps solaire vrai.
	equation_temps = 7.678 * math.sin(1.374 + (2 * pi * (jour - 81)/365)) - 9.87 * math.sin(2 * (2 * pi * (jour - 81)/365)) 
	
	# le temps solaire vrai permet de calculer la position du soleil
	temps_solaire_vrai = heure + (longitude_est/15) - fuseau - (equation_temps/60) 
	
	# l'unite de l'angle horaire doit être en radians car dans le langage python l'angle des fonctions trigonometriques est considere comme etant en radians
	angle_horaire = math.radians((temps_solaire_vrai - 12) * 15)
	
	# declinaison en radians
	declinaison = math.asin(0.4 * math.sin(math.radians(0.986 * jour - 80)))
	
	# hauteur solaire en radians
	hauteur_solaire = math.asin(math.sin(math.radians(latitude)) * math.sin(declinaison) + math.cos(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire))
	
	#Afin d'eviter des erreurs dans la suite du calcul nous limitons certaines valeurs de hauteur solaire a une valeur nulle
	if hauteur_solaire < math.radians(2):
		hauteur_solaire = 0
	
	#L'azimut est l'angle compte a partir du sud, positivement vers l'ouest, negativement vers l'est (temps_solaire_vrai < 12), entre le plan vertical du soleil a l'instant donne et le plan meridien local.
	if temps_solaire_vrai < 12:
		azimut = -1 * math.acos((math.sin(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire) - math.cos(math.radians(latitude)) * math.sin(declinaison)) / math.cos(hauteur_solaire)) 
	else:
		azimut = math.acos((math.sin(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire) - math.cos(math.radians(latitude)) * math.sin(declinaison)) / math.cos(hauteur_solaire))
	
	
	#Cette formule permet de calculer le cosinus de l'angle d'incidence
	#	(angle entre la normale au plan recepteur et le rayon solaire incident)
	# Le math.cos(0) signifie que le plan suit la trajectoire du soleil.
	cosinus_angle_incidence = math.cos(hauteur_solaire)*math.sin(0.5*pi)*math.cos(0)+math.sin(hauteur_solaire)*math.cos(0.5*pi)

	#pour eviter une division par 0 dans la formule ligne 224, lorsque la hauteur solaire est egale a 0 on considere que le flux direct normal est nul	
	if hauteur_solaire == 0:
		direct_normal = 0
	else:
		direct_normal = direct_horizontal/math.sin(hauteur_solaire)

	if diffus_horizontal == 0:
		flux_diffus_circumsolaire = 0
		flux_diffus_voute_horizon = 0	
	else:
		flux_diffus_circumsolaire = flux_diffus_circumsolaire_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence))
		# calcul de l'azimut de la face arriere de la paroi verticale en degre
		if math.degrees(azimut) < 0:
			azimut_arr = math.degrees(azimut) + 180
		else:
			azimut_arr = math.degrees(azimut) - 180
		
		azimut_arriere = math.radians(azimut_arr)
		#cosinus de l'angle d'incidence de la face arriere de la paroi verticale
		cosinus_angle_incidence_arriere = math.cos(hauteur_solaire) * math.sin(0.5 * pi) * math.cos(azimut - azimut_arriere) + math.sin(hauteur_solaire) * math.cos(0.5 * pi)
		
		flux_diffus_voute_horizon = flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence)) + flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence_arriere))
	
	if flux_diffus_circumsolaire < 0:
		flux_diffus_circumsolaire = 0

	if flux_diffus_voute_horizon < 0:
		flux_diffus_voute_horizon = 2 * diffus_horizontal * 0.5
	
	if hauteur_solaire == 0:
		flux_direct = 0
	else:
		flux_direct = direct_horizontal * cosinus_angle_incidence / math.sin(hauteur_solaire)
		
	flux_reflechi_direct = albedo * direct_horizontal * cosinus_angle_incidence	* (1 - math.cos(pi/2))/2
	flux_reflechi_diffus = albedo * diffus_horizontal * (1 - math.cos(pi/2))/2
	
	# geometrie du cylindre
	rayon_cylindre = 0.17
	hauteur_cylindre = 1.73
	# calcul des flux recus par le cylindre 
	# 	> direct = anisotrope (calcul integral)
	flux_direct_cylindre = 2 * flux_direct * rayon_cylindre * hauteur_cylindre
	# 	> diffus = anisotrope circumsolaire (calcul integral) + isotrope diffus et horizon
	flux_diffus_cylindre = (2 * flux_diffus_circumsolaire + pi * flux_diffus_voute_horizon) * rayon_cylindre * hauteur_cylindre		
	# 	> reflechi  = anisotrope (calcul integral) + isotrope
	flux_reflechi_cylindre = (2 * flux_reflechi_direct + 2 * pi * flux_reflechi_diffus) * rayon_cylindre * hauteur_cylindre	
		
	#Valeurs en Watts
	return(flux_direct_cylindre, flux_diffus_cylindre, flux_reflechi_cylindre)

def temperature_mrt_fluxsolaire(temperature_mrt_classique, diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est):
	sigma = 5.67 * 0.00000001 #constante de stefan boltzmann
	epsilon = 0.97 #emissivite du corps humain
	#temperature radiante en Kelvin
	temperature_radiante = math.pow(temperature_mrt_classique + 273.15, 4)
	f_eff = 0.75 #facteur de surface de rayonnement effectif pour un individu debout d'après JENDRITZKY and NUBLER
	alpha_clo = 0.8 #coefficient d'absorption des vêtements
	# geometrie cylindre
	rayon_cylindre = 0.17
	hauteur_cylindre = 1.73
	S_cylindre = 2*pi*rayon_cylindre*hauteur_cylindre
	
	#flux direct, diffus, reflechi avec diffus anisotrope en W/m²
	flux_direct, flux_diffus, flux_reflechi = flux_cylindre_anisotrope(diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est) / S_cylindre
	temperature_mrt_solaire = math.pow(temperature_radiante + (alpha_clo * flux_direct + f_eff * alpha_clo * flux_diffus + f_eff * alpha_clo * flux_reflechi)/(sigma * epsilon), 0.25) - 273.15
	
	return temperature_mrt_solaire
Non classé

Scilab / Windows : parallelization of the genetic algorithm optim_GA

Scilab includes a genetic algorithm in the function optim_GA.

To save time when evaluating each individual in a generation, it is possible to parallelize the evaluations: they are independent!

We propose here a modification of the Scilaboptim_ga routine allowing you to evaluate in parallel a number of individuals equal to the number of available threads on the computer (the reduction of computing time is thus proportional to the number of cores of your machine)

This operation is valid under Linux only, the Windows version of Scilab 6 not supporting parallelism: « In this current version of Scilab, parallel_run uses only one core on Windows platforms. » (source).

/////////////////////////////////////////////////////////////////////
//_______  ______    _______  _______ 
//|   _   ||    _ |  |       ||       |
//|  |_|  ||   | ||  |    ___||    _  |
//|       ||   |_||_ |   |___ |   |_| |
//|       ||    __  ||    ___||    ___|
//|   _   ||   |  | ||   |___ |   |    
//|__| |__||___|  |_||_______||___|    
// AREP - 16, av. d'Ivry, 7501/3 Paris, FRANCE
/////////////////////////////////////////////////////////////////////
// Scilab (www.scilab.org) - This file is part of Scilab
// Copyright (C) 2008 - Yann COLLETTE <yann[dot]collette[at]renault[dot]com>
// Copyright (C) 2014 - Michael Baudin <michael[dot]baudin[at]contrib[dot]scilab[dot]org>
//
//	26/6/2017 : added the "parallel_run" capability for the evaluation of individuals (Edouard Walther)
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt                                  //
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Last modification : 26/06/2017                                 //
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Contact : edouard[dot]walther[at]arep[dot]com                                
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


function [pop_opt, fobj_pop_opt, pop_init, fobj_pop_init] = optim_GA_parallel(ga_f, pop_size, nb_generation, p_mut, p_cross, Log, param)

    [nargout, nargin] = argn();

    if ~isdef("param", "local") then
        param = [];
    end

    [codage_func, err]    = get_param(param, "codage_func", coding_ga_identity);
    [init_func, err]      = get_param(param, "init_func", init_ga_default);
    [crossover_func, err] = get_param(param, "crossover_func", crossover_ga_default);
    [mutation_func, err]  = get_param(param, "mutation_func", mutation_ga_default);
    [selection_func, err] = get_param(param, "selection_func", selection_ga_elitist);
    [nb_couples, err]     = get_param(param, "nb_couples", 100);
    [pressure, err]       = get_param(param, "pressure", 0.05);
    [output_func, err] = get_param(param, "output_func", output_ga_default);

    if ~isdef("ga_f", "local") then
        error(sprintf(gettext("%s: ga_f is mandatory"), "optim_ga"));
    else
        if typeof(ga_f) == "list" then
            deff("y = _ga_f(x)", "y = ga_f(1)(x, ga_f(2:$))");
        else
            deff("y = _ga_f(x)", "y = ga_f(x)");
        end
    end

    if ~isdef("pop_size", "local") then
        pop_size = 100;
    end
    if ~isdef("nb_generation", "local") then
        nb_generation = 10;
    end
    if ~isdef("p_mut", "local") then
        p_mut = 0.1;
    end
    if ~isdef("p_cross", "local") then
        p_cross = 0.7;
    end
    if ~isdef("Log", "local") then
        Log = %F;
    end

    // Initialization of the population
    Pop = list();
    Pop = init_func(pop_size, param);

    if (nargout >= 3) then
        pop_init = Pop;
    end

    // Code the individuals
    Pop = codage_func(Pop, "code", param);

	// Getting the objective function for each individual
	//disp("Extraction list...");
	[vec_param_pop]=Pop(1);
	vec_param_total=vec_param_pop;
    for i = 2:length(Pop)
        [vec_param_pop]=Pop(i);
		vec_param_total=cat(2,vec_param_total,vec_param_pop);
    end

	// First launch 
	//disp("Execution for the first population...");
	[FObj_Pop]=parallel_run(vec_param_total,_ga_f);
	FObj_Pop=FObj_Pop';
		
    if (nargout == 4) then
        fobj_pop_init = FObj_Pop;
    end

    FObj_Pop_Max = max(FObj_Pop);
    FObj_Pop_Min = min(FObj_Pop);

    // Normalization of the efficiency
    Efficiency = (1 - pressure) * (FObj_Pop_Max - FObj_Pop) / max([FObj_Pop_Max - FObj_Pop_Min %eps]) + pressure;

	
	//disp("Starting optimisation with GA...");
	
    // The genetic algorithm
    for i = 1:nb_generation
        //
        // Selection
        //
        Indiv1 = list();
        Indiv2 = list();
        Wheel = cumsum(Efficiency);
        for j = 1:nb_couples
            // Selection of the first individual in the couple
            Shoot = grand(1, 1, "unf", 0, Wheel($));
            Index = find(Shoot <= Wheel, 1);
            Indiv1(j)      = Pop(Index);
            FObj_Indiv1(j) = FObj_Pop(Index);
            // Selection of the second individual in the couple
            Shoot = grand(1, 1, "unf", 0, Wheel($));
            Index = 1;
            Index = find(Shoot <= Wheel, 1);
            Indiv2(j)      = Pop(Index);
            FObj_Indiv2(j) = FObj_Pop(Index);
        end
        //
        // Crossover
        //
        for j = 1:nb_couples
            if (p_cross>grand(1, 1, "def")) then
                [x1, x2] = crossover_func(Indiv1(j), Indiv2(j), param);
                Indiv1(j) = x1;
                Indiv2(j) = x2;
                ToCompute_I1(j) = %T;
                ToCompute_I2(j) = %T;
            else
                ToCompute_I1(j) = %F;
                ToCompute_I2(j) = %F;
            end
        end
        //
        // Mutation
        //
        for j = 1:nb_couples
            if (p_mut>grand(1, 1, "def")) then
                x1 = mutation_func(Indiv1(j), param);
                Indiv1(j) = x1;
                ToCompute_I1(j) = %T;
            end
            if (p_mut>grand(1, 1, "def")) then
                x2 = mutation_func(Indiv2(j), param);
                Indiv2(j) = x2;
                ToCompute_I2(j) = %T;
            end
        end
        //
        // Computation of the objective functions
		
		k=0;kk=0; // counters to iterate 
		for j = 1:nb_couples // for all couples in the population
			if ToCompute_I1(j) then// if to be computed
				k=k+1;
				if k==1 then // create the first vector of parameters
					[vec_param_pop1]=Indiv1(j);
				else // concatenate for parallel_run
					[vec_param_indiv1]=Indiv1(j);
					indices_indiv1(k)=j;
					vec_param_pop1=cat(2,vec_param_pop1,vec_param_indiv1);
				end
			end
			if ToCompute_I2(j) then// if to be computed
				kk=kk+1;
				if kk==1 then
					[vec_param_pop2]=Indiv2(j);
				else
					[vec_param_indiv2]=Indiv2(j);
					indices_indiv2(kk)=j;
					vec_param_pop2=cat(2,vec_param_pop2,vec_param_indiv2);
				end
			end
		end
		
		// Parallel_run
		//disp("Parallel launch for Indiv1...");
		[objectifs_Indiv1]=parallel_run(vec_param_pop1,_ga_f);
		objectifs_Indiv1=objectifs_Indiv1';
		//disp("Parallel launch for Indiv2...");
		[objectifs_Indiv2]=parallel_run(vec_param_pop2,_ga_f);
		objectifs_Indiv2=objectifs_Indiv2';
		
		// Updating indexes
		//disp("Updating FObj1 ...");
		for k=1:length(objectifs_Indiv1)
			if indices_indiv1(k)<> 0 then
				FObj_Indiv1(indices_indiv1(k))= objectifs_Indiv1(k);
			end;
		end
		for k=1:length(objectifs_Indiv2)
			if indices_indiv2(k)<> 0 then
				FObj_Indiv2(indices_indiv2(k))= objectifs_Indiv2(k);
			end;
		end

        // Reinit ToCompute lists
        ToCompute_I1 = ToCompute_I1 & %F;
        ToCompute_I2 = ToCompute_I2 & %F;
        // Recombination
        [Pop, FObj_Pop] = selection_func(Pop, Indiv1, Indiv2, FObj_Pop, FObj_Indiv1, FObj_Indiv2, [], [], [], param);
        // Callback for plotting / printing intermediate results or stopping the algorithm
        if (Log) then
            stop = output_func(i, nb_generation, Pop, FObj_Pop, param);
            if (stop) then
                break
            end
        end
    end

    pop_opt  = Pop;
    pop_opt  = codage_func(pop_opt, "decode", param);
    fobj_pop_opt = FObj_Pop;
endfunction

 

Non classé

S.E.T. calculation in Python

Still in the comfort models, the SET* comfort index code is given below. The index is calculated according to the most up-to-date models described in the article "Comfort Modeling" on this website.

###################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
#
#  AREP, 16 avenue d'Ivry,75013 Paris, FRANCE
##################################################################
#Calcul de la SET en Python
##################################################################
# Last modification : 12/03/2018                                 #
##################################################################
# Copyright (C) 2018 Edouard Walther                             #  
# This program 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   #
# of the License, or (at your option) any later version.         #
##################################################################
# Contact : edouard[dot]walther[at]arep[dot]fr
##################################################################

import math
import random

def pv_sat(T):
	if T >= 0:
		pv_sat2 = pow(10, (2.7877 + (7.625*T)/(241 + T)))
	else :
		pv_sat2 = pow(10,(2.7877 + (9.756*T)/(272.7 + T)))
	return pv_sat2

def pv_calc(T, RH):
	pv_calc2 = (RH * pv_sat(T))/100
	return pv_calc2

def w(T, RH, p):
	w2 = 0.622 * (pv_calc(T, RH)/(p - pv_calc(T, RH)))
	return w2

def Cp_ah(T, RH, p):
	cpa = 1006
	cpv = 1830
	water = w(T, RH, p)
	Cp_ah2 = (cpa + water * cpv)/(1 + water)
	return Cp_ah2

def v_spe(T, RH, p):
	v_spe2 = (461.24 * (T + 273.15) * (0.622 + w(T, RH, p)))/p
	return v_spe2

def h_radiation(T_rad, T_surf):
	h_radiation2 = 0.72 * 5.67 * 0.00000001 * ((T_rad + T_surf) + 2 * 273.15) * (pow((T_rad + 273.15), 2) + pow((T_surf + 273.15), 2)) * 0.97 #<--- emissivite = 0,97
	return h_radiation2


###############################################################################
############################ FONCTION SET* ####################################
############################################################################### 	

def modele_metabolique_SET(RadTempMtx, WindSpeedMtx, T_ambient, phi_ambient, p_ambient, hauteur, masse, fat, Cst_dilat,Cst_sweat, Cst_constr, T_core_set, T_skin_set, SkinBloodFlow_set,U_muscle_fat_skin, C_shiv):
	T_skin = T_skin_set
	T_core = T_core_set
	dT = 60
	#p_ambient = 101325
	#met = 1.1
	#clo = 1
	#i_m = 0.45
	# metabolisme masculin en W
	age=30
	surface =  0.203*pow(hauteur,0.725)*pow(masse,0.425) #Surface exterieure du sujet [m2] 
	#genre=random.random()
	genre = 0.1 # homme
	if genre<0.5:
		metabolisme_W = 3.45 * math.pow(masse, 0.75) * (1.0 + 0.004 * (30.0 - age) + 0.01 * (hauteur * 100.0 / math.pow(hauteur, 1.0 / 3.0) - 43.4))
	else:
		metabolisme_W = 3.19 * math.pow(masse, 0.75) * (1.0 + 0.004 * (30.0 - age) + 0.018* (hauteur * 100.0 / math.pow(hauteur, 1.0 / 3.0) - 42.1))
	met = metabolisme_W/surface/58.2
	#print surface, met
	SkinBloodFlow = 6.3
	minutes_metab = 60
	minutes = minutes_metab - 1	
	duree=(minutes + 1) * 60
	temps=0
	while temps < duree: 
		#variables fonctionnelles 
		compteur = 0 #variable de boucle
		temps = temps + dT
		i_m = 0.45
		i_m_static = i_m
		clo = 1 #0.155 m2.K/W
		clo_static = clo
		#evolution metabolique dynamique
		#if temps < duree/2:
		#	met = 2.2
		#	v_walk = 4/3.6
		#else:
		#	met = 1.2
		#	v_walk = 0  	
		  	#WindSpeedMtx = 0.2	
		  	#clo dynamique
		#if v_walk < 0.7:
		#	v_walk = 0.0052 * (met * 58.2 - 58)
		#corr_T = math.exp(0.042 - 0.398 * WindSpeedMtx + 0.066 * WindSpeedMtx**2 - 0.378 * v_walk + 0.094 * v_walk**2)	
		#if WindSpeedMtx > 3.5:
		#	corr_T = 0.582
		# ATTENTION : ICI CHANGER "WindSpeedMtx" car sinon on le change a chaque pas de temps !
		#WindSpeedMtx = math.sqrt(WindSpeedMtx**2 + v_walk**2)
		#clo = clo_static * corr_T
		#i_m = i_m_static * (4.9 - 6.5 * corr_T + 2.6 * corr_T**2)
		#Constantes du corps humain
		KCLO = 0.25  #coefficient augmentation surface d'echange
#		masse = vect[1] #masse moyenne sur une population
#		R_muscle_fat_skin = 5.28
		#Cp_body = 0.97*3600 #capacite calorifique du corps humain [J/(kg.K)]
		body_mass = masse
#		fat = 15 # pourcentage masse graisseuse
		fat_mass = fat/100*body_mass
		Cp_body = fat_mass/body_mass*2510 + (body_mass - fat_mass)/body_mass*3650 # Modele HAVENITH
	  	#Constantes de regulation de l'organisme
		SBC = 0.0000000567 #Constante de Stefann-Boltzmann
#		Cst_sweat = 170
#		Cst_dilat = 200
#		Cst_constr = 0.5
		#Valeurs de consignes de la regulation du corps humain
#		T_skin_set = 33.7 #temperature de peau
#		T_core_set = 36.8 #temperature interne
		T_body_set = 0.1*T_skin_set + 0.9*T_core_set #temperature corporelle
#		SkinBloodFlow_set = 6.3 #debit sanguin [L/m2.h]
		#Conversion d'unites
		p_atmosphere = p_ambient/1000 #conversion en kPa
		p_atmosphere = p_atmosphere*0.009869 #conversion en atm
		R_clo = 0.155*clo
		#correction de la veture  ########## 
		if clo < 0.5:
			f_surf_clo = 1 + 0.2*clo
			
		else:
			f_surf_clo = 1 + 0.15*clo
		
		#calcul du nombre de Lewis
		Lewis = 2434 * v_spe(T_ambient, phi_ambient*100, p_ambient)/(Cp_ah(T_ambient, phi_ambient*100, p_ambient) * 1.04 * pow(0.83, (2/3))) * (18/8.32/(T_ambient + 273.15))
		#Calculs initiaux du metabolisme
		RM = met * 58.2
		Metab = met * 58.2
		w_crit = 0.59 * pow(WindSpeedMtx, (-0.08))
		#Calcul des coefficients d'echange convectif
		h_c = 3 * pow(p_atmosphere, 0.53)
		h_c_vent = 8.600001 * pow(WindSpeedMtx * p_atmosphere, 0.53)
		h_c = max(h_c, h_c_vent)
		#Coefficient d'echange radiatif
		h_r = 4.7
		#Coefficient d'echange global
		h_g = h_r + h_c
		#Resistance thermique convective+radiative
		R_air = 1/(f_surf_clo * h_g)
		#Temperature operative
		T_op = (h_r * RadTempMtx + h_c * T_ambient)/h_g #Sous forme de matrice		
		#Temperature superficielle de veture
		T_clo = T_op + (T_skin - T_op)/(h_g * (R_air + R_clo))
		T_clo_OLD = T_clo + 0.5	      
	############### Boucle calcul T_clo ###########################################
		while abs(T_clo - T_clo_OLD) > 0.001:
	  		T_clo_OLD = T_clo
	  		h_r = h_radiation(RadTempMtx, T_clo)# Sous forme de matrice
	  		h_g = h_r + h_c
	  		R_air = 1/(f_surf_clo * h_g)
	  		T_op = (h_r * RadTempMtx + h_c * T_ambient)/h_g# Sous forme de matrice
	  		T_clo = (R_air * T_skin + R_clo * T_op)/(R_air + R_clo)
	  		compteur = compteur + 1
	  		if compteur > 20:
	  			break
	###############################################################################
	  	#SkinBloodFlow = SkinBloodFlow_set
	  	#Temperature corporelle
	  	#alpha = 0.0417737 + 0.7451833/(SkinBloodFlow + 0.585417)
	  	#T_body = alpha * T_skin + (1 - alpha) * T_core	  			
	################ Modele de regulation du corps humain #######################
		#Skin signal
		signal_skin = T_skin - T_skin_set
		if signal_skin > 0:
			warm_skin = signal_skin
			cold_skin = 0
		else :
	  		warm_skin = 0
	  		cold_skin = -signal_skin
		# Core signal
		signal_core = T_core - T_core_set
		if signal_core > 0:
			warm_core = signal_core
			cold_core = 0
		else:
			warm_core = 0
			cold_core = -signal_core
		#Debit sanguin
		SkinBloodFlow = (SkinBloodFlow_set + Cst_dilat * warm_core)/(1 + Cst_constr * cold_skin)
		if SkinBloodFlow > 90 :
			SkinBloodFlow = 90
		if SkinBloodFlow < 0.5 :
			SkinBloodFlow = 0.5
		#Temperature corporelle
		alpha = 0.0417737 + 0.7451833/(SkinBloodFlow + 0.585417)
		T_body = alpha * T_skin + (1 - alpha) * T_core
		#Corps/Body
		signal_body = T_body - T_body_set
		if signal_body > 0 :
			warm_body = signal_body
			cold_body = 0
		else :
			warm_body = 0
			cold_body = -signal_body
		#Debit sudation
		qm_sweat = Cst_sweat * warm_body * math.exp(warm_skin/10.7)
		if qm_sweat > 500:
			qm_sweat = 500
		#Chaleur latente maximale echangee
		R_vap_tot = (R_clo + R_air)/(Lewis * i_m)
		E_max = (pv_sat(T_skin) - phi_ambient * pv_sat(T_ambient))/R_vap_tot
		h_e = (2.2 * h_c)/(1 + 0.928 * R_clo * h_c)/133.322
		#E_max = h_e * (pv_sat(T_skin) - phi_ambient * pv_sat(T_ambient))
		#Chaleur latente sudation
		E_sweat = 0.68 * qm_sweat
		#Ratios
		pcent_sweat = E_sweat/E_max #Part d'energie echangee due a la sudation adimensionnel
		pcent_wet = 0.06 + 0.94 * pcent_sweat #Part de la surface du corps mouille
		#Chaleur latente echangee par diffusion de l'eau a travers la couche cutanee
		E_diff = pcent_wet * E_max - E_sweat
		#Chaleur latente totale echangee par la peau
		E_skin = E_sweat + E_diff
		#Sudation supercritique
		if pcent_wet > w_crit :
	  		pcent_wet = w_crit
	  		pcent_sweat = w_crit/0.94
	  		E_sweat = pcent_sweat * E_max
	  		E_diff = 0.06 * (1 - pcent_sweat) * E_max
	  		E_skin = E_sweat + E_diff
	  		drip_cond_nope = 1
		#Condensation (la pression de vapeur saturante a la temperature de la peau est inferieure a la pression de vapeur de l'air ambiant
		elif E_max < 0 :
	  		E_diff = 0
	  		E_sweat = 0
	  		pcent_wet = w_crit
	  		pcent_sweat = w_crit
	  		E_skin = E_max
	  		drip_cond_nope = -2
		else:
			drip_cond_nope = 0
		w_skin = pcent_wet
		#Frisson/shivering
		M_shiv = C_shiv * cold_skin * cold_core
		Metab = RM + M_shiv
		w_sweat_global = E_sweat/E_skin #Part d'humidite due a la sudation
		#Flux echange entre l'interieur du corps (noyau) et la peau
		Flx_core_skin = (T_core - T_skin) * (U_muscle_fat_skin + 1.163 * SkinBloodFlow) # 1.163 = 4 200/3 600 * 1000
		#Chaleur sensible echangee entre la peau et l'exterieur
		DRY = (T_skin - T_op)/(R_air + R_clo)
		#Chaleur sensible echangee par la respiration
		T_expir = 32.5 + 0.066 * T_ambient + 1.99 * 0.000001 * phi_ambient * pv_sat(T_ambient)
		C_resp = 0.0014 * Metab * (T_expir - T_ambient) #
		#Chaleur latente echangee par la respiration
		E_resp = 0.000017251 * Metab * (pv_sat(35.5) - phi_ambient * pv_sat(T_ambient))
		
		#Accumulation d'energie par la peau
		SSK = Flx_core_skin - DRY - E_skin
		Accumulation_skin = SSK
		
		#Accumulation d'energie par le corps
		SCR = Metab - Flx_core_skin - E_resp - C_resp
		Accumulation_core = SCR
		
		#Modification des temperatures par l'effet de l'accumulation
		
		#Methode 1
		dT_skin = SSK * surface * dT /(alpha * masse * Cp_body)
		dT_core = SCR * surface * dT / ((1-alpha) * masse * Cp_body)
		
		#Methode 2
		TCSK = Cp_body/3600 * alpha * masse
		TCCR = Cp_body/3600 * (1 - alpha) * masse
		
		DTSK = (SSK * surface)/(TCSK * dT)
		DTCR = (SCR * surface)/(TCCR * dT)
	
		T_skin = T_skin + DTSK
		T_core = T_core + DTCR
		T_body = alpha * T_skin + (1 - alpha) * T_core		
		#Energie totale echangee par la peau
		H_skin = DRY + E_skin
		#Metabolisme net
		RN = Metab
		#Chaleur latente echangee par sudation en etat de confort
		E_comf = 0.42 * (RN - (1*58.2))
		if E_comf < 0:
			E_comf = 0
		E_max = E_max * w_crit
		#Chaleur latente evaporative requise pour la thermoregulation
		E_req = RN - E_resp - C_resp - DRY
		E_sweat_global = E_sweat
		#Coefficient d'echange chaleur sensible
		HD = 1/(R_air + R_clo)
		#Coefficient d'echange evaporatif
		HE = 1/R_vap_tot
		#Pression de vapeur saturante a la temperature cutanee
		PSSK = pv_sat(T_skin)
		#Coefficients d'echange
		h_r_SET = h_r
		if met < 0.85:
			h_c_SET = 3
		else:
			h_c_SET = 5.66 * pow((met - 0.85), 0.39)
			if h_c_SET < 3:
				h_c_SET = 3
		h_g_SET = h_c_SET + h_r_SET
		#CLO metabolique
		RCLOS = 1.52/(met + 0.6944) - 0.1835
		#Resistance thermique de la veture
		RCLS = 0.155 * RCLOS
		#Correction de la surface d'echange due a la veture
		f_surf_clo_SET = 1 + KCLO * RCLOS
		#Facteur d'efficacite de BURTON
		F_clo_SET = 1/(1 + 0.155 * f_surf_clo_SET * h_g_SET * RCLOS)
		#Index de permeabilite de la veture
		i_m_SET = 0.38
		#Indice de permeation de la veture
		i_clo_SET = i_m_SET * h_c_SET/h_g_SET * (1 - F_clo_SET)/(h_c_SET/h_g_SET - F_clo_SET * i_m_SET)
		#Resistance convective + radiative corrigee
		R_air_SET = 1/(f_surf_clo_SET * h_g_SET)
		#Resistances evaporatives
		REAS = 1/(Lewis * f_surf_clo_SET * h_c_SET)
		RECLS = RCLS/(Lewis * i_clo_SET)
		#Resistance totale au transfert de chaleur sensible
		HD_S = 1/(R_air_SET + RCLS)
		HE_S = 1/(REAS + RECLS)
		#Variables de resolution de SET/ET
		Delta = 0.0001
		dx = 100
		X_OLD = T_skin - H_skin/HD_S
		while abs(dx) > 0.001:
			#compteurr=compteurr+1
			ERR1 = H_skin - HD_S * (T_skin - X_OLD) - w_skin * HE_S * (PSSK - 0.5 * pv_sat(X_OLD))
			ERR2 = H_skin - HD_S * (T_skin - (X_OLD + Delta)) - w_skin * HE_S * (PSSK - 0.5 * pv_sat(X_OLD + Delta))
			#if ERR2==ERR1:
			#	print("attention sortie brutale de boucle...")
			#	break
			x = X_OLD - Delta * ERR1/(ERR2 - ERR1)
			dx = x - X_OLD
			X_OLD = x
		ET_global = x
	return ET_global
Non classé

The corrected P.E.T.

Here is the modified P.E.T. comfort index, as described in thefollowing article. The modifications are taken into account.

##################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
# AREP - 16, av. d'Ivry, 75013 Paris, FRANCE
##################################################################
# based on: Peter Hoeppe PET fortran code, from:
# "Urban climatic map and standards for wind environment - Feasibility study, Technical Input Report No.1",
# Edouard Walther and Quentin Goestchel
# Most of the changes in the implementaion are explained in the resolution function comments                                  #
##################################################################
# Last modification : 10/04/2018                                 #
##################################################################
# Copyright (C) 2018 Édouard WALTHER                             #  
# This program 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   #
# of the License, or (at your option) any later version.         #
##################################################################
# Contact : edouard[dot]walther[at]arep[dot]com                                
##################################################################


import os
import numpy as np
import math as math
import scipy.optimize as optimize

## Implementation of the skin and core temperatures set values #######
tc_set=36.6 # 36.8
tsk_set=34 # 33.7
tbody_set=0.1*tsk_set+0.9*tc_set # Calculation of the body temperature through a weighted average


## Skin blood flow calculation function: #######
def vasoC(tcore,tsk):
    #Set value signals to consider every cases:
    sig_skin = tsk_set - tsk
    sig_core = tcore - tc_set
    if sig_core<0:
        # In this case, Tcore<Tc_set: the body needs to keep the heat --> the blood flow is reduced
        sig_core=0.
    if sig_skin<0:
        # In this case, Tsk>Tsk_set: the body needs to loose heat --> the blood flow is increased
        sig_skin=0.
    # 6.3 L/m^2/h is the set value of the blood flow
    qmblood = (6.3 + 75. * sig_core) / (1. + 0.5 * sig_skin)
    # 90 L/m^2/h is the blood flow upper limit, not sustainable for a human being
    if qmblood>90:
        qmblood=90.
    # Alpha can be used to calculate tbody in ameliorated models
    #alfa = 0.04177 + 0.74518 / (qmblood + 0.585417)
    alfa = 0.1
    return (qmblood,alfa)


## Sweating flow calculation function: #######
def Suda(tbody,tsk):
    sig_body = tbody - tbody_set
    sig_skin = tsk - tsk_set
    if sig_body<0:
        #In this case, Tbody<Tbody_set: the body needs to keep the heat --> The sweat flow is 0
        sig_body=0.
    if sig_skin<0:
        # In this case, Tsk<Tsk_set: the body needs to keep the heat --> the sweat flow is reduced
        sig_skin=0.
    #qmsw = 170 * sig_body * math.exp((sig_skin) / 10.7)  # [g/m2/h] Expression from Gagge Model
    qmsw = 304.94*10**(-3) * sig_body
    # 90 L/m^2/h is the blood flow upper limit, not sustainable for a human being
    if qmsw > 500:
        qmsw = 500
    return (qmsw)


## Vectorial MEMI balance calculation function for the 3 node model: #######
def Syst(T, Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,mode):
    ## Conversion of T vector in an array
    arr = np.ones((3,1))
    arr[0,0]=T[0] #Corresponds to T_core
    arr[1,0]=T[1] #Corresponds to T_skin
    arr[2,0]=T[2] #Corresponds to T_clothes
    T=arr
    enbal_vec = np.zeros((3,1)) #Useful for the vectorial expression of the balance

    ## Area parameters of the body: #######
    Adu = 0.203 * mbody ** 0.425 * ht ** 0.725
    feff=0.725
    if pos == 1 or pos == 3:
        feff = 0.725
    if pos == 2:
        feff = 0.696
    # Calculation of the Burton coefficient, k = 0.31 for Hoeppe:
    fcl = 1 + (0.31 * icl) # Increase of the heat exchange surface  depending on clothing level
    facl = (173.51 * icl - 2.36 - 100.76 * icl * icl + 19.28 * icl ** 3.0) / 100
    Aclo = Adu * facl + Adu * (fcl - 1.0)
    Aeffr = Adu * feff  # Effective radiative area depending on the position of the subject
    # Partial pressure of water in the air depending on relative humidity and air temperature:
    if mode: # actual environment
        vpa = HR / 100.0 * 6.105 * math.exp(17.27 * Ta / (237.7 + Ta )) #[hPa]
    else:# mode==False implies we are calculating the PET
        vpa= 12 # [hPa] vapour pressure of the standard environment

    # Convection coefficient depending on wind velocity and subject position
    hc = 0
    if pos == 1:
        hc = 2.67 + (6.5 *v**0.67)
    if pos == 2:
        hc = 2.26 + (7.42 *v**0.67)
    if pos == 3:
        hc = 8.6 * (v ** 0.513)
    hc = hc * (p / po) ** 0.55

    # Basic metabolism for men and women in [W] #######
    metab_female = 3.19 * mbody**0.75 * (1.0 + 0.004 * (30.0 - age) + 0.018 * (ht * 100.0 / mbody**(1.0/3.0) - 42.1))
    metab_male = 3.45 * mbody**0.75 * (1.0 + 0.004 * (30.0 - age) + 0.01 * (ht * 100.0 / mbody**(1.0/3.0) - 43.4))
    # Source term : metabolic activity
    eswpot = (M + metab_male)/Adu
    fec = (M + metab_female)/Adu
    he = 0.0
    # Attribution of internal energy depending on the sex of the subject
    if sex == 1:
        he = eswpot
    elif sex == 2:
        he = fec
    h = he *(1.0 - eta)  # [W/m2]

    # Respiratory energy losses
    # Expired air temperature calculation:
    texp = 0.47 * Ta + 21.0  # [degC]

    # Pulmonary flow rate
    dventpulm = he * 1.44 * 10.0**(-6.0)

    # Sensible heat energy loss:
    eres = cair * (Ta - texp) * dventpulm  # [W/m2]

    # Latent heat energy loss:
    vpexp = 6.11 * 10.0**(7.45 * texp / (235.0 + texp))
    erel = 0.623 * Lvap / p * (vpa-vpexp) * dventpulm  # [W/m2]
    ere = eres + erel  # [W/m2]

    # Clothed fraction of the body approximation
    rcl = icl / 6.45  # Conversion in m2.K/W
    y = 0
    if facl > 1.0:
        facl = 1.0
    if icl >= 2.0:
        y = 1.0
    if icl > 0.6 and icl < 2.0:
        y = (ht - 0.2)/ht
    if icl <= 0.6 and icl > 0.3:
        y = 0.5
    if icl <= 0.3 and icl > 0.0:
        y = 0.1
    # calculation of the closing radius depending on the clothing level (6.28 = 2* pi !)
    r2 = Adu * (fcl - 1.0 + facl) / (6.28 * ht * y)  # External radius
    r1 = facl * Adu /(6.28 * ht * y)  # Internal radius
    di = r2 - r1

    # Calculation of the thermal resistance of the body
    alpha = vasoC(T[0,0],T[1,0])[1]
    tbody = alpha * T[1,0] + (1 - alpha) * T[0,0]
    htcl = (6.28 * ht * y * di)/(rcl * math.log(r2 / r1)*Aclo)  # [W/(m2.K)]

    # Calculation of sweat losses

    qmsw = Suda(tbody,T[1,0])
    # Lvap Latent heat of evaporation : 2400 [J/g] divided by 3600 for [g/m2/h] to [g/m2/s]
    esw = 2400 * qmsw / 3600  # [W/m2]
    # Saturation vapor pressure at temperature Tsk and for 100% HR
    Pvsk = 6.105 * math.exp((17.27 * (T[1,0]+273.15) - 4717.03)/ (237.7 + T[1,0])) # [hPa]

    rscl=0.155*icl
    Lw = 16.7 * 10 ** (-1)  # [K/hPa] Lewis factor
    he_diff = hc * Lw
    fecl=1/(1+0.92*hc*rscl)
    emax = he_diff * fecl * (Pvsk - vpa)
    w = esw / emax  # Dermal wetness
    if w > 1:
        w=1
        delta = esw-emax
        if delta < 0:
            esw=emax
    if esw < 0:
        esw=0
    i_m=0.3
    R_ecl=(1/(fcl*hc) + rscl)/(Lw*i_m)
    #R_ecl=0.79*1e7 #version Hoeppe
    ediff = (1 - w)*(Pvsk - vpa)/R_ecl  # Basal perspiration
    evap = -(ediff + esw)  # [W/m2]

    # Radiation losses
    # For bare skin area:
    rbare = Aeffr*(1.0 - facl) * emsk * sigm * ((Tmrt + 273.15)**(4.0) - (T[1,0] + 273.15)**(4.0))/Adu
    # For dressed area:
    rclo = feff * Aclo * emcl * sigm * ((Tmrt + 273.15)**(4.0) - (T[2,0] + 273.15)**(4.0))/Adu
    rsum = rclo+rbare


    ## Convection losses #######
    cbare = hc * (Ta - T[1,0]) * Adu * (1.0 - facl)/Adu  # [w/m^2]
    cclo = hc * (Ta - T[2,0]) * Aclo/Adu  # [W/m^2]
    csum = cclo+cbare

    ## Balance equations of the 3-nodes model
    enbal_vec[0,0] = h + ere - (vasoC(T[0,0],T[1,0])[0]/3600*cb+5.28)*(T[0,0]-T[1,0]) # Core balance [W/m^2]
    enbal_vec[1,0] = rbare + cbare + evap + (vasoC(T[0,0],T[1,0])[0]/3600*cb+5.28)*(T[0,0]-T[1,0]) - htcl*(T[1,0]-T[2,0])  # Skin balance [W/m^2]
    enbal_vec[2,0] = cclo + rclo + htcl*(T[1,0]-T[2,0]) # Clothes balance [W/m^2]
    enbal_scal = h + ere + rsum + csum +evap

    if mode:
        return [enbal_vec[0,0],enbal_vec[1,0],enbal_vec[2,0]] #List of the balances values
    else:
        return enbal_scal #Scalar balance used in the PET calculation

# Solving the system
def resolution(Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl, Tx):
    Tn = optimize.fsolve(Syst,Tx ,args=(Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,True))
    return (Tn, 1)


# PET calculation with dichotomy method 
def PET (age, sex, ht, mbody, pos, M, icl, Tstable,a,b,eps):
    # Definition of a function with the input variables of the PET reference situation
    def f(Tx):
        return Syst(Tstable, Tx, Tx, 50, 0.1, age, sex, ht, mbody, pos, M, 0.9,False)
    Ti = a # Start index of the browsing interval
    Tf = b # End index of the browsing interval
    pet = 0
    while Tf-Ti>eps: # Dichotomy loop
        if f(Ti)*f(pet)<0:
            Tf = pet
        else:
            Ti = pet
        pet = (Ti + Tf) / 2.0
    return pet

# Input data
# Environment constants
po = 1013.25 #[hPa]
rob = 1.06 # Blood density kg/L
cb = 3.64 * 1000. # Blood specific heat [j/kg/k]
cair = 1.01 * 1000. # Air specific heat  [J./kg/K-]
emsk = 0.99 # Skin emissivity
emcl = 0.95 # Clothes emissivity
Lvap = 2.42 * 10. ** 6. # Latent heat of evaporation [J/Kg]
sigm = 5.67 * 10. ** (-8.) # Stefan-Boltzmann constant [W/(m2*K^(-4))]
eta = 0. # Body efficiency

# Initialisation of Temperature vector with respectively: Tcore, Tskin, Tcl
T = [38,40,40]
eps = 10**(-3) #Numerical precision

# Dichotomy browsning parameters
a = -40
b = 60
# Input data of the PET 
Ta=50 # Air temperature in oC
Tmrt=50 # Mean radiant temperature in oC
HR=50 # Air relative humidity %
v=0.1 # Wind velocity m/s
age = 35
sex = 1 # 1 for men and 2 for women
pos = 1
mbody = 75 #[Kg]
ht = 1.80 #[m]
p = 1013.25 #[hPa]
M_activity = 80 #[W] Activity
icl = 0.9# [clo] Clothing level

#initialisation pour le premier calcul
Tstable = resolution(Ta,Tmrt,HR,v,age,sex,ht,mbody,pos,M_activity,icl,T)[0]
print(PET(age, sex, ht, mbody, pos, M_activity, icl, Tstable, -30, 90, eps))
Non classé

Running Energyplus in parallel on Windows

A very useful shell script for parallel executions ofEnergyPlus on a Windows system, with the command-line Git-Bash (!)

#!/bin/bash
##################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
#
##################################################################
# ce script en bash lance energypus en parallele 
#(tout est dans le "&" esperluette)                              #
##################################################################
# Last modification : 01/04/2018                                 #
##################################################################
# Copyright (C) 2018 Alexis SAUVAGEON                            #  
# This program 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   #
# of the License, or (at your option) any later version.         #
##################################################################
# Contact : alexis[dot]sauvageon[at]arep[dot]com                                
##################################################################


cd ${0%/*} || exit 1 # un peu de magie au debut (tire de openfoam)
# define nbr cpus
ncpus=$1
ncpus=30

if [ -z "${1}" ]; then
    	ncpus=30
fi
# pour tout fichier "f" dans le dossier courant "*/"
for f in Trajectoire*/ 
do # faire :
	echo '----- lancement energyplus pour ' $f \n # ecrire qu'on commence avec le fichier "f"
	# lancer energyplus (meme cmd que dans python)
	g=${f::-1} # j'enleve le "/" final pour eviter sa redondance dans la ligne suivante avec celui de "/Simulation_*.idf"
	# on execute E+ pour tous les dossiers au meme niveau que le script :
	#		> IDF de chaque dossier 
	#		> EPW au meme niveau que le script
	exec C:/EnergyPlusV8-3-0/energyplus.exe -d $f -w ./in.epw $g/Simulation_*.idf &
	#		> on check si E+ tourne et on l'affiche dans le terminal
	process=$(ps cax | grep energyplus)
	if [ $? -eq 0 ]; then
		echo "Process is running."
	else
		echo "Process is not running."
	fi
	#		> on compte le nombre de processus actif E+
	numproc=$(ps -ef | grep -v grep | grep energyplus | wc -l)
	#		> tant que nbr E+ actif > ncpu autorise, on attend et on met a jour numproc
  while [[ $numproc -ge $ncpus ]]
  do
    sleep 5
    numproc=$(ps -ef | grep -v grep | grep energyplus | wc -l)
    echo $numproc 
  done
done