Trimarchi, Daniele

English

Southampton (e-Prints Soton)

OpenFOAM WorkshopSouthampton, 21 October 2011         UNIVERSITY OFSouthamptonSchool of Engineering SciencesProgramming session: from the C++ basics to the compilation of user librariesDaniele Trimarchidaniele.trimarchi@soton.ac.ukPresentation overview:• C++ basics- main, declarations, types- memory management: references and pointers- Object orientation: ‣ classes‣ inheritance- header and source ﬁles• Compiling applications - make, wmake, wmake libso. Linking libraries.• Writing simple applications in OpenFOAM- general structure (includes, skim between time-dirs...)- examples: divergence, INDT30 m15 mPresentation overview:• Modifying existing libraries- overview of the force class- modifying the class‣ Change the source and modify the names‣ compile the library‣ call the library during the execution• Adding classes- FSInterface class- Linking classes trough pointers25 m20 mPART 1 C++ basicsC++ basicsHello world code... it only prints a message on screen C++ basicsDeclare variables, assign values and output the resultsC++ basicsDeﬁne and use of the functionsC++ basicsReferences and pointersiA rB rCComputer memory: sequential representation& iA & rB & rC• References ( & ) are the addresses of the variables in the computer memory.  References are constant. • The value of a reference (ex: 0x7fff5fbfb604) can be stored in a particular type of variable, called pointer.• The  pointer  is  a  variable  of  the  same  type  of  the variable it points to. Integer pointer points to integer variable!C++ basicsint * APtr   =  & iAﬂoat * rPtr =  & rB...rPtr = & rCDeclaration AssignationRe-Assignation• Using references and pointers  is  convenient: it avoids the need of making copies when passing arguments to functions (for large arrays this is very convenient!)ﬂoat rB = 1.1;ﬂoat rC = 0.6; ResAdd = Add(rB, rC);cout << ResAdd << endl;mainﬂoat c = a + b;functioncopy: a = rBcopy: b = rCcopy: ResAdd = crB rC resAdd a b cMemory allocated for main Memory allocated for functionC++ basics• The deferencing operator * restitutes the value of the variable pointed by the pointer. So: ﬂoat rB = 1.1;ﬂoat * rPtr =  & rBcout<<   rB <<endl;           1.1cout<< * rPtr <<endl;        1.1cout<< & rB<<endl;          0x7fff5fbfb604Variable declaration and assignation:Pointer declaration and assignation: Similar syntax, but very different meaning!!Pointer deferencing: C++ basicsOptimizing the functionPass references to the function(no copies!)Assign the pointer:ﬂoat * a = & rB12Deference and return therefore  a ﬂoat value3Assign the returned value to the ﬂoat variable 4C++ basics• CRectangle classClass declarationPrivate membersPublic membersand functions Public functiondeﬁnitionDeclarationsDeﬁnitionObject declarationCall the function member of the classC++ basicsThe function is called with the operator  “ . ”  :  rect.setValues(3,4)The function can also be called using a pointer: CRectangle tria(2,3);CRectangle * triaPtr;triaPtr = & tria;cout << tria.area() << endl;cout << * tria.area() << endl;cout << triaPtr -≻ area() << endldeclare the variabledeclare the pointerassign the pointercall fcn using objectcall fcn using object, by de-referencing pointercall fcn using pointerC++ basics• The constructor is a class-member function called when the object is initially build• Dual of the Constructor is the Destructor• As for every function, different arguments deﬁne different constructors; the constructor can also call other functions Standard syntax for the constructorC++ basics• Inheritance: Cpolygon classCPolygonCRectangle CTrianglewidth, heightArea AreaInherited class: class CRectangle : public CPolygon... is a ... Every object has its own function deﬁnition. The function behaves differently accordingly to the type of the object !!C++ basics• Main (.cpp; .C), Header (.H) and Source (.C) ﬁlesDeclarations.H ﬁleDeﬁnition.C ﬁleMain code.cpp, .C ﬁleC++ basics• Main (.cpp; .C), Header (.H) and Source (.C) ﬁlesDeclarations.H ﬁleDeﬁnition.C ﬁleMain code.cpp, .C ﬁleC++ basics• Compiling applications: passing from human readable instructions to binaries• Unix environment: call makeﬁle trough “make”Set the compilerCompiler optionsCompiler Flags; external librariesSource ﬁlesExecutable• The “make” command is overwritten in OpenFOAM by “wmake” and “wmake libso”ReferencesSeveral books are available on C++. For example: Deitel, DeitelC++, How to programed. Prentice HallA book is generally better (expecially for C++)...but on-line useful guidance can be also found:http://www.cplusplus.com/doc/tutorial/PART 2 Writing simple applications in OpenFOAMGeneral structure of an OpenFOAM application...A lot of nice code... src/ﬁniteVolume/cfdTools/general/include/fvCFD.Hseveral other include: classes for time, mesh, geometry, math constants...src/OpenFOAM/include/createTime.HDeclares runTime, object of the class Foam::Time. Constructor deﬁned in src/OpenFOAM/db/Time/Time.H, line193General structure of an OpenFOAM application...A lot of nice code...  File in the source directoryDeclares a VolScalarField called divergence to be written in every time-step folder Calculating the divergence ∇·∇ ·￿ v =∂v1x1+∂v2x2+∂v3x3= δikvi,kCheck for existing time dirsOpenFOAM version of the for loop. Equivalent to: for(int timeI=0; i<timeDirs.size(); timeI++){...}Declares and reads the ﬁeld U from the selected time directorycalculates the ﬁeldCompiling the applicationSolutionDivergence.CcreateFields.HMakeﬁlesoptionsWorking directory:solutionDivergence.CEXE = $(FOAM_USER_APPBIN)/DivUEXE_INC = \    -I$(LIB_SRC)/ﬁniteVolume/lnIncludeEXE_LIBS = -lﬁniteVolumeSource ﬁleCompile applicationinclude headers for FVinclude FV libraryThe application is compiled typing at terminal the command: wmakeCalculating the Normalised invariant of the deformation tensorD =SijSij − WijWijSijSij + WijWij;Sij =12∂ui∂xj+∂uj∂xi;Wij =12∂ui∂xj−∂uj∂xiCreateFields.h Core of the main codeTextDouble inner product operator, see Programmer guide P23PART 3Modifying existing librariesConstructor of the class forcesPath:OpenFOAM/src/postProcessing/functionObjects/forces/forces/Force.CReference to the Object Registry. This is a list of the entities pertaining to an objectReference to the controlDictCall the member function forces::readRead the entries in the controlDictReading entries from the controlDictsystem/controlDict:functions(  forces  {    type forces;    functionObjectLibs ("libforces.dylib");                                                                                                                                                                 outputControl outputTime;    patches (wing);                                                                                                                                                                                               pName p;    Uname U;    rhoName rhoInf;    rhoInf 1.2; //Reference density                                                                                                                                                                                                                                    pRef 0;    CofR (0 0 0); //Origin for moments                                                                                                                                                                                                                                               })patches on which forces will be integrated Calculating the forcesThe virtual function write() is called during the execution. This calls forces::calcForcesMoment(), where the calculation is performedOpenFOAM iterator. it corresponds to a for cyclemesh is object of the class fvMeshThe expression returns a vector with the cell centres of the chosen patch Calculating the forcesThe virtual function write() is called during the execution. This calls forces::calcForcesMoment(), where the calculation is performedSfb is the (reference to) the face area vectorIt is here multiplied for the pressure boundaryField => pf returns the vector of forces on the chosen patchRe-compiling the forces libraryThe basic idea of the openFOAM environement is: find something similar and modify it as you like, but DO NOT TOUCH THE ORIGINAL SOURCES!STEP 1: copy the forces directory from the original location into another directorySTEP 2: copy also the Make folderSTEP 3: substitute strings and modify (all) the file names( sed ‘s/forces/Myforces/g’ forces.C > Myforces.C )STEP 4: modify the local functionObject.H file (add the new class to the list of loadable functions )#include "Myforces.H"   ...namespace Foam   {     typedef OutputFilterFunctionObject<Myforces>forcesFunctionObject;   }   ...Re-compiling the forces librarySTEP 4: Modify the Make/files:STEP 5: modify the Make/options file:include all what was needed by the original library!STEP 6: compile with wmake libsoEXE_INC = \.... all what was already there ...    -I$(LIB_SRC)/postProcessing/functionObjects/forces/lnIncludeMyforces.CforcesFunctionObject.CLIB = $(FOAM_USER_LIBBIN)/LibMyforcesRe-compiling the forces librarySTEP 7: Add the entries in the controlDict, in order for the library to be loaded and used:libs ( "libMyforces.dylib" ) ;                  Load the library (.dylib on MAC; .so on Linux)     functions (                                                                   Use the library   Myforces            Search in the library                         {                          for the entry called     type Myforces;     functionObjectLibs ("libMyforces.dylib");     outputControl outputTime;      patches (BottWall);       pName p;     Uname U;     rhoName rhoInf;     rhoInf 1.0;      pRef 0;      CofR (0 0 0);    });ReferencesSeveral examples can be found on-line:http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2007/http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2008/http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2009/ http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2010/ See for example the work of A. Gonzales.But have a deep look at the whole web-site, there’s a lot of enlightening material!!PART 4Adding new classes in OpenFOAMthe FSInterface classScope of this class• The class is designed to manage the mesh motions. It is used in the main of the solver pimpleDyMFOAM:• establish the communication (if needed) with the external solver trough MPI• send pressure data and retrieve mesh displacement data• communicate with the AitkenControl class, in charge for calculating the dynamic relaxation factor: • move the ﬂuid mesh (ALE framework)uk+1 = ωk ˜ uk+1 +( 1− ωk)ukScope of this classPimpleDyMFOAM External solver FSInterface MPIC++...whatever works......C++, fortran...DataOut (p)DataIn (u)AitkenControlAitPtrincludeincludeMultiple Program Multiple Data type environment, the external solver is “spawned” during the execution time. This generates a communicator we can use for exchanging data (white arrows)Use of the classIn the main solver: Include and declare new classes deﬁnitionFoam classes needed for the mesh motionDeclare the new objects using the right arguments!Constructor of FSInterface:header (.H) ﬁlesource (.C) ﬁle... a lot of other stuff...References are constant: they MUST be initialised at the creation. In this case this is done by passing the values to the constructorArguments passed to the classAssign values to the class membersExecute other fcnslinking to the AitkenControl class:An object of the type AitkenControl is created right before the object FSInterface.  A reference to this object is passed to the constructor. This reference is stored in a pointermain code:alfa2 is instantiated in the constructor, the pointer is referenced also in the constructor:Searching the mesh motion entriesDepending upon the motion solver, the mesh motion is stored in a ﬁeld “pointDisplacements”, “cellDisplacement” or “motionU”. A ﬁeld in OpenFOAM is deﬁned as: internalField + boundaryFieldImposing the motion of a boundary means writing the motion in the  correspondent  entry  of  the  boundaryField.  For  example, “&MeshDisplacement” is the address of the BoundaryField, while “&mDisp” is the address of the mesh interface we want to move Searching the mesh motion entriesIOdictionary couplingDict(IOobject    (	" CouplingDict",	 runTime.constant(),	 mesh,	 IOobject::MUST_READ,	 IOobject::NO_WRITE        ));word temp(couplingDict.lookup("ﬂuidPatch"))word interface = temp;label ﬂuidPatchID = mesh.boundaryMesh().ﬁndPatchID(interface);Opens the FOAM dictionary named “couplingDict”, to be read only. The dict must be placed in the folder constantSearches and reads the entry “ﬂuidPatch” in CouplingDictSearches the the entry “ﬂuidPatchName” in the mesh boundaryField. Returns a (integer) label: the id of the patch in the mesh orderSearching the mesh motion entriesWe need now to ﬁnd the references to the mesh motion boundaryField. Using displacementLaplacian...pointVectorField & PointDisplacement = const_cast<pointVectorField&>(mesh.objectRegistry::lookupObject<pointVectorField>("pointDisplacement"));pDispPtr = & ( refCast<vectorField>(PointDisplacement.boundaryField()[ﬂuidPatchID]));Search in the objectRegistry of the mesh an object of the type: pointVectorField called pointDisplacement. Whatever its type, impose (const_cast) to be a reference of the type PointVectorFieldThe entry number “ﬂuidPatchID” is the reference to the interface mesh motion. Store this reference into a pointerPointDisplacement.boundaryField()[ ﬂuidSideI ] == U_kp1The mesh motion is imposed using the surcharged operator == :  Searching the mesh motion entriesWe need now to ﬁnd the references to the mesh motion boundaryField. Using laplaceFaceDecomposition...const fvMesh& motionMesh =runTime.objectRegistry::lookupObject<fvMesh>(movingRegionName);tetPointVectorField& motionU =const_cast<tetPointVectorField&>( motionMesh.objectRegistry::lookupObject<tetPointVectorField>("motionU") );	tetPolyMesh& tetMesh = const_cast<tetPolyMesh&>(motionU.mesh());motionUFluidPatchPtr = &refCast<ﬁxedValueTetPolyPatchVectorField>(  motionU.boundaryField()[ﬂuidPatchID]  );   * motionUFluidPatchPtr == ( U_kp1 - U_old ) / runTime.deltaT().value();The mesh motion is imposed using the surcharged operator == :  Scheme of the classAlgorithmFLUID STRUCTUREdt = dt +1pkukdt = dtImplicit  coupling:  the  equilibrium  within  the  time-step  is veriﬁed  using  a  ﬁxed-point  between  the  ﬂuid  and  the structural solver. Although this algorithm is computationally expensive, it is unconditionally stableThis is realised in OpenFOAM using a while loop within the time-step.  The convergence criterion is veriﬁed by a function of the class FSInterfaceAlgorithm (main)Aitken iterations looptime loopLaunch the structural solverSend the number of iteration to the interface (different calculations if iterCntr == 1)Send the pressures, retreive the displacementsConvergence is checked by interfaceCompilingAs for every FOAM application, we need to edit: FSIﬂuidFOAMtet.CEXE = $(FOAM_USER_APPBIN)/FSIﬂuidFOAMtetMPIMake/ﬁles:Source ﬁleName and path of the compiled applicationMake/options:include $(RULES)/mplib$(WM_MPLIB)EXE_INC = \...all remains as in the original ﬁle...EXE_LIBS = \...all remains as in the original ﬁle, but add:    -lmpi \    $(WM_DECOMP_LIBS)And ﬁnally some results!FLUID DENSITYFLUID VISCOSITYSTRUCTURE’sTHICKNESSSTRUCTURE’s DENSITYSTRUCTURE’sPOISSON RATTIOSTRUCTURE’sYOUNG’S MODULUSρf νf ts ρs νs EsKg/m3 m2/s m Kg/m3 - N/m21.0 0.01 0.002 500 0 250u =1− cos(2πt/5)

Programming session: from the C++ basics to the compilation of user libraries

Abstract

Similar works

Full text

Available Versions

Southampton (e-Prints Soton)