exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Merge branch 'geometricImmersedBoundary' into development

Browse source
This commit is contained in:
Hrvoje Jasak 2018-02-07 14:40:02 +00:00
commit bb8423f3e2
520 changed files with 34237 additions and 9273 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
.gitattributes vendored Normal file
View file

@ -0,0 +1,3 @@
*.H gitlab-language=cpp
*.C gitlab-language=cpp

38
ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch vendored
View file

@ -1,6 +1,31 @@
diff -ruN ParMGridGen-1.0_orig/Makefile ParMGridGen-1.0/Makefile
--- ParMGridGen-1.0_orig/Makefile 2017-04-04 15:02:44.020713666 +0200
+++ ParMGridGen-1.0/Makefile 2017-04-04 15:21:48.582647336 +0200
@@ -1,16 +1,21 @@
default:
+ (mkdir bin)
(cd MGridGen ; make)
serial:
+ (mkdir bin)
(cd MGridGen ; make)
parallel:
+ (mkdir bin)
(cd MGridGen ; make)
(cd ParMGridGen ; make)
clean:
+ (mkdir bin)
(cd MGridGen ; make clean)
(cd ParMGridGen ; make clean )
realclean:
+ (mkdir bin)
(cd MGridGen ; make realclean )
(cd ParMGridGen ; make realclean )
diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
--- ParMGridGen-1.0_orig/Makefile.in 2001-12-04 16:30:33.000000000 -0800
+++ ParMGridGen-1.0/Makefile.in 2013-08-22 20:07:33.491171127 -0700
--- ParMGridGen-1.0_orig/Makefile.in 2017-04-04 15:02:44.012713543 +0200
+++ ParMGridGen-1.0/Makefile.in 2017-04-04 15:06:00.159742074 +0200
@@ -1,6 +1,6 @@
#--------------------------------------------------------------------------
# Which make to use
@ -18,6 +43,15 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
# Which loader to use
LD = cc
@@ -22,7 +22,7 @@
LDOPTIONS = -O3
# Where to put the executable
-BINDIR = ../..
+BINDIR = ../../bin
# Additional libraries
DMALLOCDIR = /usr/local
@@ -33,22 +33,25 @@
# In which directories to look for any additional libraries

38
ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch_darwin vendored
View file

@ -79,8 +79,8 @@ diff -ruN ParMGridGen-1.0_orig/MGridGen/Programs/Makefile ParMGridGen-1.0/MGridG
ifeq ($(ddmalloc),yes)
DEBUGFLAGS := $(DEBUGFLAGS) -DDMALLOC -DDEBUG
diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
--- ParMGridGen-1.0_orig/Makefile.in 2011-12-24 13:54:44.000000000 -0500
+++ ParMGridGen-1.0/Makefile.in 2011-12-24 13:49:26.000000000 -0500
--- ParMGridGen-1.0_orig/Makefile.in 2001-12-05 01:30:33.000000000 +0100
+++ ParMGridGen-1.0/Makefile.in 2017-04-04 15:36:04.695980033 +0200
@@ -1,6 +1,6 @@
#--------------------------------------------------------------------------
# Which make to use
@ -98,6 +98,15 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
# Which loader to use
LD = cc
@@ -22,7 +22,7 @@
LDOPTIONS = -O3
# Where to put the executable
-BINDIR = ../..
+BINDIR = ../../bin
# Additional libraries
DMALLOCDIR = /usr/local
@@ -33,22 +33,25 @@
# In which directories to look for any additional libraries
@ -129,3 +138,28 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
#--------------------------------------------------------------------------
#
diff -ruN ParMGridGen-1.0_orig/Makefile ParMGridGen-1.0/Makefile
--- ParMGridGen-1.0_orig/Makefile 2001-11-09 00:41:22.000000000 +0100
+++ ParMGridGen-1.0/Makefile 2017-04-04 14:51:04.033914737 +0200
@@ -1,16 +1,21 @@
default:
+ (mkdir bin)
(cd MGridGen ; make)
serial:
+ (mkdir bin)
(cd MGridGen ; make)
parallel:
+ (mkdir bin)
(cd MGridGen ; make)
(cd ParMGridGen ; make)
clean:
+ (mkdir bin)
(cd MGridGen ; make clean)
(cd ParMGridGen ; make clean )
realclean:
+ (mkdir bin)
(cd MGridGen ; make realclean )
(cd ParMGridGen ; make realclean )

12
ThirdParty/rpmBuild/SOURCES/ParaView-4.4.0.patch_darwin vendored Normal file
View file

@ -0,0 +1,12 @@
diff -ruN ParaView-4.3.1_orig/Applications/ParaView-4.3.1_extra_install_Darwin.cmake ParaView-4.3.1/Applications/ParaView-4.3.1_extra_install_Darwin.cmake
--- ParaView-4.3.1_orig/Applications/ParaView-4.3.1_extra_install_Darwin.cmake 1969-12-31 19:00:00.000000000 -0500
+++ ParaView-4.3.1/Applications/ParaView-4.3.1_extra_install_Darwin.cmake 2013-10-02 19:00:00.000000000 -0400
@@ -0,0 +1,8 @@
+#
+# Additional install rules for Mac OS X platforms
+#
+INSTALL (DIRECTORY buildObj/bin/paraview.app
+ DESTINATION ${CMAKE_INSTALL_PREFIX}/bin
+ USE_SOURCE_PERMISSIONS
+ COMPONENT Runtime)
+

16
ThirdParty/rpmBuild/SOURCES/libccmio-2.6.1.patch_0 vendored
View file

@ -1,13 +1,13 @@
diff -ruN libccmio-2.6.1_orig/config/config.gnu.to.star libccmio-2.6.1/config/config.gnu.to.star
--- libccmio-2.6.1_orig/config/config.gnu.to.star 2007-12-17 18:32:03.000000000 -0500
+++ libccmio-2.6.1/config/config.gnu.to.star 2010-12-18 14:50:35.000000000 -0500
--- libccmio-2.6.1_orig/config/config.gnu.to.star 2007-12-18 00:32:03.000000000 +0100
+++ libccmio-2.6.1/config/config.gnu.to.star 2017-04-04 13:41:35.000000000 +0200
@@ -1,4 +1,4 @@
-#!/bin/sh
+#!/bin/bash
# $Id: config.gnu.to.star,v 1.4 2006/06/05 21:12:16 geoffp Exp $
@@ -34,6 +34,12 @@
@@ -34,6 +34,13 @@
x86_64-unknown-linux-gnu-null) echo linux64_2.4-x86-glibc_2.2.5 ;;
ppc64-unknown-linux-gnu-null) echo linux64_2.6-pwr4-glibc_2.3.3 ;;
i386-apple-darwin8-null) echo i386-apple-darwin8 ;;
@ -17,6 +17,7 @@ diff -ruN libccmio-2.6.1_orig/config/config.gnu.to.star libccmio-2.6.1/config/co
+ i386-apple-darwin13-null) echo i386-apple-darwin13 ;;
+ i386-apple-darwin14-null) echo i386-apple-darwin14 ;;
+ i386-apple-darwin15-null) echo i386-apple-darwin15 ;;
+ i386-apple-darwin16-null) echo i386-apple-darwin16 ;;
*) echo unknown ;;
esac
@ -30,15 +31,15 @@ diff -ruN libccmio-2.6.1_orig/config/config.guess libccmio-2.6.1/config/config.g
# Copyright (C) 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001
# Free Software Foundation, Inc.
diff -ruN libccmio-2.6.1_orig/config/config.system libccmio-2.6.1/config/config.system
--- libccmio-2.6.1_orig/config/config.system 2008-02-25 22:07:16.000000000 -0500
+++ libccmio-2.6.1/config/config.system 2010-12-18 14:51:34.000000000 -0500
--- libccmio-2.6.1_orig/config/config.system 2008-02-26 04:07:16.000000000 +0100
+++ libccmio-2.6.1/config/config.system 2017-04-04 13:42:15.000000000 +0200
@@ -1,4 +1,4 @@
-#! /bin/sh
+#! /bin/bash
# $Id: config.system,v 1.2 2005/09/29 22:19:19 geoffp Exp $
@@ -87,6 +87,27 @@
@@ -87,6 +87,30 @@
i386-apple-darwin8.11.1)
echo i386-apple-darwin8 ;;
@ -62,6 +63,9 @@ diff -ruN libccmio-2.6.1_orig/config/config.system libccmio-2.6.1/config/config.
+
+ i386-apple-darwin15.* )
+ echo i386-apple-darwin15 ;;
+
+ i386-apple-darwin16.* )
+ echo i386-apple-darwin16 ;;
+
*)
echo unknown

1343
ThirdParty/rpmBuild/SOURCES/mesquite-2.1.2_patch1 vendored Normal file
View file

File diff suppressed because it is too large Load diff

1
ThirdParty/rpmBuild/SPECS/libccmio-2.6.1.spec vendored
View file

@ -105,6 +105,7 @@ Patch0: libccmio-2.6.1.patch_0
[ ! -d config/i386-apple-darwin13 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin13
[ ! -d config/i386-apple-darwin14 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin14
[ ! -d config/i386-apple-darwin15 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin15
[ ! -d config/i386-apple-darwin16 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin16
%endif
# Warning:
# 1: The name of the ADF library will be renamed to libadf_ccmio since this

2
ThirdParty/rpmBuild/SPECS/mesquite-2.1.2.spec vendored
View file

@ -80,6 +80,7 @@ Source: %url/%{name}-%{version}.tar.gz
Prefix: %{_prefix}
Group: Development/Tools
Patch0: mesquite-2.1.2_patch0
Patch1: mesquite-2.1.2_patch1
%define _installPrefix %{_prefix}/packages/%{name}-%{version}/platforms/%{_WM_OPTIONS}
@ -90,6 +91,7 @@ Patch0: mesquite-2.1.2_patch0
%setup -q
%patch0 -p1
%patch1 -p1
%build
# export WM settings in a form that GNU configure recognizes

5
applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C
View file

@ -72,6 +72,11 @@ int main(int argc, char *argv[])
U.internalField() = vector::zero;
}
p.correctBoundaryConditions();
U.correctBoundaryConditions();
phi == (fvc::interpolate(U) & mesh.Sf());
# include "volContinuity.H"
# include "meshCourantNo.H"

16
applications/solvers/basic/sixDOFSolver/sixDOFSolver.C
View file

@ -36,7 +36,7 @@ Description
#include "objectRegistry.H"
#include "foamTime.H"
#include "ODESolver.H"
#include "sixDOFbodies.H"
#include "sixDOFBodies.H"
#include "OFstream.H"
using namespace Foam;
@ -48,9 +48,12 @@ int main(int argc, char *argv[])
# include "setRootCase.H"
# include "createTime.H"
sixDOFbodies structure(runTime);
sixDOFBodies structure(runTime);
OFstream of(runTime.path()/"motion.dat");
// Write header for output file
of << "# Time, CoG_x, CoG_y, CoG_z, omega_x, omega_y, omega_z" << endl;
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
@ -73,11 +76,14 @@ int main(int argc, char *argv[])
<< structure()[bodyI].omegaAverage().value() << nl
<< "Current omega in time instant = "
<< structure()[bodyI].omega().value() << nl
<< "Average omegaAbs in time step = "
<< structure()[bodyI].omegaAverageAbsolute().value() << nl
<< endl;
of << structure()[bodyI].X().value().x() << tab;
of << structure()[bodyI].X().value().x() << tab
<< structure()[bodyI].X().value().y() << tab
<< structure()[bodyI].X().value().z() << tab
<< structure()[bodyI].omega().value().x() << tab
<< structure()[bodyI].omega().value().y() << tab
<< structure()[bodyI].omega().value().z() << tab;
}
of << endl;

49
applications/solvers/immersedBoundary/SOLVER_COOKBOOK
View file

@ -1,49 +0,0 @@
Make/options: add
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I../immersedBoundary/lnInclude
-ltriSurface \
-lmeshTools \
-L$(FOAM_USER_LIBBIN) -limmersedBoundary \
createFields.H: add
# include "createIbMasks.H"
solver:
1) add immersed boundary headers
#include "immersedBoundaryFvPatch.H"
#include "immersedBoundaryAdjustPhi.H"
2) on calculation of face fluxes (or their components), mask with faceIbMask
3) before to adjustPhi, add:
// Adjust immersed boundary fluxes
immersedBoundaryAdjustPhi(phi, U);
4) on explicit updates, add correctBoundaryConditions();
eg.
U = fvc::reconstruct(phi);
U.correctBoundaryConditions();
5) On reports of continuity error, add masking:
Info<< "IB-masked continuity error = "
<< mag(cellIbMask*fvc::div(phi))().weightedAverage(mesh.V()).value()
<< endl;
or use immersedBoundaryContinuityErrs.H
5) Chenge Courant number check to be IB-sensitive, using
immersedBoundaryCourantNo.H

3
applications/solvers/immersedBoundary/icoDyMIbFoam/Make/files
View file

@ -1,3 +0,0 @@
icoDyMIbFoam.C
EXE = $(FOAM_APPBIN)/icoDyMIbFoam

56
applications/solvers/immersedBoundary/icoDyMIbFoam/createFields.H
View file

@ -1,56 +0,0 @@
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar nu
(
transportProperties.lookup("nu")
);
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, pimple.dict(), pRefCell, pRefValue);
mesh.schemesDict().setFluxRequired(p.name());

152
applications/solvers/immersedBoundary/icoDyMIbFoam/icoDyMIbFoam.C
View file

@ -1,152 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
icoDyMOversetFoam
Description
Transient solver for incompressible, laminar flow of Newtonian fluids
with dynamic mesh and immersed boundary mesh support.
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "immersedBoundaryFvPatch.H"
#include "immersedBoundaryAdjustPhi.H"
#include "pimpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createDynamicFvMesh.H"
pimpleControl pimple(mesh);
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
// Make the fluxes absolute
fvc::makeAbsolute(phi, U);
bool meshChanged = mesh.update();
reduce(meshChanged, orOp<bool>());
Info<< "Mesh update" << meshChanged << endl;
# include "createIbMasks.H"
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
# include "CourantNo.H"
// Pressure-velocity corrector
while (pimple.loop())
{
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}
// --- PISO loop
while (pimple.correct())
{
volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();
// Immersed boundary update
U.correctBoundaryConditions();
phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf());
// Adjust immersed boundary fluxes
immersedBoundaryAdjustPhi(phi, U);
adjustPhi(phi, U, p);
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve
(
mesh.solutionDict().solver
(
p.select(pimple.finalInnerIter())
)
);
if (pimple.finalNonOrthogonalIter())
{
phi -= pEqn.flux();
}
}
# include "immersedBoundaryContinuityErrs.H"
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

3
applications/solvers/immersedBoundary/icoIbFoam/Make/files
View file

@ -1,3 +0,0 @@
icoIbFoam.C
EXE = $(FOAM_APPBIN)/icoIbFoam

13
applications/solvers/immersedBoundary/icoIbFoam/Make/options
View file

@ -1,13 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lsurfMesh \
-lsampling \
-ldynamicMesh \
-limmersedBoundary \
-llduSolvers

139
applications/solvers/immersedBoundary/icoIbFoam/icoIbFoam.C
View file

@ -1,139 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
icoIbFoam
Description
Transient solver for incompressible, laminar flow of Newtonian fluids
with immersed boundary support.
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "immersedBoundaryFvPatch.H"
#include "immersedBoundaryAdjustPhi.H"
#include "pimpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
pimpleControl pimple(mesh);
# include "createIbMasks.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;
# include "CourantNo.H"
// Pressure-velocity corrector
while (pimple.loop())
{
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}
// --- PISO loop
while (pimple.correct())
{
volScalarField rUA = 1.0/UEqn.A();
U = rUA*UEqn.H();
// Immersed boundary update
U.correctBoundaryConditions();
phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf());
// Adjust immersed boundary fluxes
immersedBoundaryAdjustPhi(phi, U);
adjustPhi(phi, U, p);
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve
(
mesh.solutionDict().solver
(
p.select(pimple.finalInnerIter())
)
);
if (pimple.finalNonOrthogonalIter())
{
phi -= pEqn.flux();
}
}
# include "immersedBoundaryContinuityErrs.H"
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

3
applications/solvers/immersedBoundary/interIbFoam/Make/files
View file

@ -1,3 +0,0 @@
interIbFoam.C
EXE = $(FOAM_APPBIN)/interIbFoam

27
applications/solvers/immersedBoundary/interIbFoam/UEqn.H
View file

@ -1,27 +0,0 @@
surfaceScalarField muEff
(
"muEff",
twoPhaseProperties.muf()
+ fvc::interpolate(rho*turbulence->nut())
);
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::laplacian(muEff, U)
- (fvc::grad(U) & fvc::grad(muEff))
);
UEqn.relax();
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
- gh*fvc::grad(rho)
- fvc::grad(pd)
);
}

56
applications/solvers/immersedBoundary/interIbFoam/alphaEqn.H
View file

@ -1,56 +0,0 @@
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
// New formulation of phir: Co condition
surfaceScalarField phir
(
"phir",
faceIbMask*interface.cAlpha()*interface.nHatf()*
min
(
scalar(0.5)/ // This is ref compression Co number for cAlpha = 1
(
mesh.surfaceInterpolation::deltaCoeffs()*
mesh.time().deltaT()
),
mag(phi/mesh.magSf())
+ dimensionedScalar("small", dimVelocity, 1e-3)
)
);
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1)
+ fvm::div(phi, alpha1, alphaScheme)
+ fvm::div
(
-fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
alpha1,
alpharScheme
)
);
alpha1Eqn.relax();
alpha1Eqn.solve();
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(alpha1) = " << min(cellIbMask*alpha1).value()
<< " Max(alpha1) = " << max(cellIbMask*alpha1).value()
<< endl;
rhoPhi = faceIbMask*(alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2);
// Limit alpha to handle IB cells
// alpha1.max(0);
// alpha1.min(1);
// Update of interface and density
interface.correct();
twoPhaseProperties.correct();
// Calculate density using limited alpha1
rho = twoPhaseProperties.rho();
rho.correctBoundaryConditions();
}

58
applications/solvers/immersedBoundary/interIbFoam/correctPhi.H
View file

@ -1,58 +0,0 @@
{
# include "continuityErrs.H"
wordList pcorrTypes
(
pd.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
for (label i=0; i<pd.boundaryField().size(); i++)
{
if (pd.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", pd.dimensions(), 0.0),
pcorrTypes
);
dimensionedScalar rUAf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
// Adjust immersed boundary fluxes
immersedBoundaryAdjustPhi(phi, U);
adjustPhi(phi, U, pcorr);
mesh.schemesDict().setFluxRequired(pcorr.name());
while(pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rUAf, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pdRefCell, pdRefValue);
pcorrEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();
}
}
# include "immersedBoundaryContinuityErrs.H"
}

118
applications/solvers/immersedBoundary/interIbFoam/createFields.H
View file

@ -1,118 +0,0 @@
Info<< "Reading field pd\n" << endl;
volScalarField pd
(
IOobject
(
"pd",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field alpha1\n" << endl;
volScalarField alpha1
(
IOobject
(
"alpha1",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
Info<< "Reading transportProperties\n" << endl;
twoPhaseMixture twoPhaseProperties(U, phi, "alpha1");
const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
// Need to store rho for ddt(rho, U)
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT
),
alpha1*rho1 + (scalar(1) - alpha1)*rho2,
alpha1.boundaryField().types()
);
rho.oldTime();
// Mass flux
// Initialisation does not matter because rhoPhi is reset after the
// alpha1 solution before it is used in the U equation.
surfaceScalarField rhoPhi
(
IOobject
(
"rho*phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
rho1*phi
);
Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("gh", g & mesh.Cf());
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pd + rho*gh
);
label pdRefCell = 0;
scalar pdRefValue = 0.0;
setRefCell(pd, pimple.dict(), pdRefCell, pdRefValue);
mesh.schemesDict().setFluxRequired(pd.name());
mesh.schemesDict().setFluxRequired(alpha1.name());
// Construct interface from alpha1 distribution
interfaceProperties interface(alpha1, U, twoPhaseProperties);
// Construct incompressible turbulence model
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
);

24
applications/solvers/immersedBoundary/interIbFoam/limitU.H
View file

@ -1,24 +0,0 @@
{
scalar limitMagU = readScalar(pimple.dict().lookup("limitMagU"));
volScalarField magU(mag(U));
scalar maxMagU = max(magU).value();
Info<< "mag(U): max: " << maxMagU
<< " avg: " << magU.weightedAverage(mesh.V()).value();
if (maxMagU > limitMagU)
{
U.internalField() *=
neg(magU.internalField() - limitMagU)
+ pos(magU.internalField() - limitMagU)*
limitMagU/(magU.internalField() + SMALL);
U.correctBoundaryConditions();
Info << " ...limiting" << endl;
}
else
{
Info << endl;
}
}

53
applications/solvers/immersedBoundary/interIbFoam/pEqn.H
View file

@ -1,53 +0,0 @@
{
volScalarField rUA = 1.0/UEqn.A();
surfaceScalarField rUAf = fvc::interpolate(rUA);
U = rUA*UEqn.H();
// Immersed boundary update
U.correctBoundaryConditions();
# include "limitU.H"
surfaceScalarField phiU
(
"phiU",
faceIbMask*(fvc::interpolate(U) & mesh.Sf())
);
phi = phiU
+ faceIbMask*
(
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
- ghf*fvc::snGrad(rho)
)*rUAf*mesh.magSf();
// Adjust immersed boundary fluxes
immersedBoundaryAdjustPhi(phi, U);
adjustPhi(phi, U, pd);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pdEqn
(
fvm::laplacian(rUAf, pd, "laplacian(rAU,pd)") == fvc::div(phi)
);
pdEqn.setReference(pdRefCell, pdRefValue);
pdEqn.solve
(
mesh.solutionDict().solver(pd.select(pimple.finalInnerIter()))
);
if (pimple.finalNonOrthogonalIter())
{
phi -= pdEqn.flux();
}
}
// Explicitly relax pressure
pd.relax();
U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
U.correctBoundaryConditions();
}

3
applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/files Normal file
View file

@ -0,0 +1,3 @@
pimpleDyMIbFoam.C
EXE = $(FOAM_APPBIN)/pimpleDyMIbFoam

27
applications/solvers/immersedBoundary/interIbFoam/Make/options → applications/solvers/immersedBoundary/pimpleDyMIbFoam/Make/options
View file

@ -1,22 +1,23 @@
EXE_INC = \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
EXE_LIBS = \
-linterfaceProperties \
-lfiniteVolume \
-limmersedBoundary \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-ldynamicMesh \
-lmeshTools \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-lmeshTools \
-lsurfMesh \
-lsampling \
-ldynamicMesh \
-limmersedBoundary \
-llduSolvers
-llduSolvers \
-L$(MESQUITE_LIB_DIR) -lmesquite

26
applications/solvers/immersedBoundary/pimpleDyMIbFoam/UEqn.H Normal file
View file

@ -0,0 +1,26 @@
// Convection-diffusion matrix
fvVectorMatrix HUEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff()
);
// Time derivative matrix
fvVectorMatrix ddtUEqn(fvm::ddt(U));
// Get under-relaxation factor
scalar UUrf =
mesh.solutionDict().equationRelaxationFactor(U.select(pimple.finalIter()));
if (pimple.momentumPredictor())
{
// Solve momentum predictor
solve
(
ddtUEqn
+ relax(HUEqn, UUrf)
==
- fvc::grad(p),
mesh.solutionDict().solver((U.select(pimple.finalIter())))
);
}

31
applications/solvers/immersedBoundary/pimpleDyMIbFoam/correctPhi.H Normal file
View file

@ -0,0 +1,31 @@
{
volScalarField pcorr("pcorr", p);
// Initialise flux with interpolated velocity
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, pcorr);
mesh.schemesDict().setFluxRequired(pcorr.name());
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(1/aU, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);
pcorrEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();
}
// Fluxes are corrected to absolute velocity and further corrected
// later. HJ, 6/Feb/2009
}
# include "continuityErrs.H"
}

11
applications/solvers/immersedBoundary/pimpleDyMIbFoam/createControls.H Normal file
View file

@ -0,0 +1,11 @@
#include "createTimeControls.H"
bool correctPhi
(
pimple.dict().lookupOrDefault("correctPhi", false)
);
bool checkMeshCourantNo
(
pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
);

43
applications/solvers/immersedBoundary/icoIbFoam/createFields.H → applications/solvers/immersedBoundary/pimpleDyMIbFoam/createFields.H
View file

@ -1,23 +1,3 @@
Info<< "Reading transportProperties\n" << endl;
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar nu
(
transportProperties.lookup("nu")
);
Info<< "Reading field p\n" << endl;
volScalarField p
(
@ -54,3 +34,26 @@
scalar pRefValue = 0.0;
setRefCell(p, pimple.dict(), pRefCell, pRefValue);
mesh.schemesDict().setFluxRequired(p.name());
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
Info<< "Reading field aU if present\n" << endl;
volScalarField aU
(
IOobject
(
"aU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
1/runTime.deltaT(),
zeroGradientFvPatchScalarField::typeName
);

58
applications/solvers/immersedBoundary/pimpleDyMIbFoam/pEqn.H Normal file
View file

@ -0,0 +1,58 @@
{
p.boundaryField().updateCoeffs();
// Prepare clean Ap without time derivative contribution and
// without contribution from under-relaxation
// HJ, 26/Oct/2015
aU = HUEqn.A();
// Store velocity under-relaxation point before using U for the flux
// precursor
U.storePrevIter();
U = HUEqn.H()/aU;
phi = (fvc::interpolate(U) & mesh.Sf());
// Global flux balance
adjustPhi(phi, U, p);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(1/aU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve
(
mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
);
if (pimple.finalNonOrthogonalIter())
{
phi -= pEqn.flux();
}
}
// Explicitly relax pressure for momentum corrector except for last corrector
if (!pimple.finalIter())
{
p.relax();
}
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
# include "movingMeshContinuityErrs.H"
U = UUrf*
(
1.0/(aU + ddtUEqn.A())*
(
U*aU - fvc::grad(p) + ddtUEqn.H()
)
)
+ (1 - UUrf)*U.prevIter();
U.correctBoundaryConditions();
}

78
applications/solvers/immersedBoundary/interIbFoam/interIbFoam.C → applications/solvers/immersedBoundary/pimpleDyMIbFoam/pimpleDyMIbFoam.C
View file

@ -22,19 +22,17 @@ License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
interIbFoam
pimpleDyMFoam.C
Description
Solver for 2 incompressible, isothermal immiscible fluids using a VOF
(volume of fluid) phase-fraction based interface capturing approach,
with immersed boundary support
The momentum and other fluid properties are of the "mixture" and a single
momentum equation is solved.
Transient solver for incompressible, flow of Newtonian fluids
with dynamic mesh using the PIMPLE (merged PISO-SIMPLE) algorithm.
Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
For a two-fluid approach see twoPhaseEulerFoam.
Consistent formulation without time-step and relaxation dependence by Jasak
Support for immersed boundary
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved
@ -42,32 +40,30 @@ Author
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "interfaceProperties.H"
#include "twoPhaseMixture.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "dynamicFvMesh.H"
#include "pimpleControl.H"
#include "immersedBoundaryPolyPatch.H"
#include "immersedBoundaryFvPatch.H"
#include "immersedBoundaryAdjustPhi.H"
#include "emptyFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createDynamicFvMesh.H"
pimpleControl pimple(mesh);
# include "readGravitationalAcceleration.H"
# include "initContinuityErrs.H"
# include "createFields.H"
# include "createIbMasks.H"
# include "createTimeControls.H"
# include "correctPhi.H"
# include "CourantNo.H"
# include "setInitialDeltaT.H"
# include "createFields.H"
# include "createControls.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -75,19 +71,46 @@ int main(int argc, char *argv[])
while (runTime.run())
{
# include "readTimeControls.H"
# include "immersedBoundaryCourantNo.H"
# include "readControls.H"
# include "CourantNo.H"
# include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// Pressure-velocity corrector
bool meshChanged = mesh.update();
U.correctBoundaryConditions();
p.correctBoundaryConditions();
aU.correctBoundaryConditions();
phi.correctBoundaryConditions();
turbulence->correct();
# include "updateIbMasks.H"
# include "volContinuity.H"
if (runTime.outputTime())
{
volScalarField divMeshPhi("divMeshPhi", mag(fvc::surfaceIntegrate(mesh.phi())));
divMeshPhi.write();
}
if (checkMeshCourantNo)
{
# include "meshCourantNo.H"
}
// Fluxes will be corrected to absolute velocity
// HJ, 6/Feb/2009
# include "correctPhi.H"
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);
// --- PIMPLE loop
while (pimple.loop())
{
twoPhaseProperties.correct();
# include "UEqn.H"
// --- PISO loop
@ -96,15 +119,6 @@ int main(int argc, char *argv[])
# include "pEqn.H"
}
# include "immersedBoundaryContinuityErrs.H"
# include "limitU.H"
p = pd + rho*gh;
p.correctBoundaryConditions();
# include "alphaEqn.H"
turbulence->correct();
}

5
applications/solvers/immersedBoundary/pimpleDyMIbFoam/readControls.H Normal file
View file

@ -0,0 +1,5 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault("correctPhi", false);
checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false);

3
applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/files
View file

@ -1,3 +0,0 @@
porousSimpleIbFoam.C
EXE = $(FOAM_APPBIN)/porousSimpleIbFoam

20
applications/solvers/immersedBoundary/porousSimpleIbFoam/Make/options
View file

@ -1,20 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools \
-lsurfMesh \
-lsampling \
-ldynamicMesh \
-limmersedBoundary \
-limmersedBoundaryTurbulence \
-llduSolvers

12
applications/solvers/immersedBoundary/porousSimpleIbFoam/UEqn.H
View file

@ -1,12 +0,0 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff()
);
UEqn().relax();
pZones.addResistance(UEqn());
solve(UEqn() == -fvc::grad(p));

44
applications/solvers/immersedBoundary/porousSimpleIbFoam/createFields.H
View file

@ -1,44 +0,0 @@
Info << "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info << "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
mesh.schemesDict().setFluxRequired(p.name());
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, laminarTransport)
);
porousZones pZones(mesh);

41
applications/solvers/immersedBoundary/porousSimpleIbFoam/pEqn.H
View file

@ -1,41 +0,0 @@
{
volScalarField AU = UEqn().A();
U = UEqn().H()/AU;
// Immersed boundary update
U.correctBoundaryConditions();
UEqn.clear();
phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf());
// Adjust immersed boundary fluxes
immersedBoundaryAdjustPhi(phi, U);
adjustPhi(phi, U, p);
// Non-orthogonal pressure corrector loop
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(1.0/AU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi -= pEqn.flux();
}
}
# include "immersedBoundaryContinuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
// Momentum corrector
U -= fvc::grad(p)/AU;
U.correctBoundaryConditions();
}

3
applications/solvers/immersedBoundary/potentialIbFoam/Make/files
View file

@ -1,3 +0,0 @@
potentialIbFoam.C
EXE = $(FOAM_APPBIN)/potentialIbFoam

16
applications/solvers/immersedBoundary/potentialIbFoam/Make/options
View file

@ -1,16 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lsurfMesh \
-lsampling \
-ldynamicMesh \
-ldynamicFvMesh \
-limmersedBoundary \
-llduSolvers

50
applications/solvers/immersedBoundary/potentialIbFoam/createFields.H
View file

@ -1,50 +0,0 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
);
// Do not reset p
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Do not reset U
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::interpolate(U) & mesh.Sf()
);
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
mesh.schemesDict().setFluxRequired(p.name());

231
applications/solvers/immersedBoundary/potentialIbFoam/potentialIbFoam.C
View file

@ -1,231 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
potentialIbFoam
Description
Potential flow solver with immersed boundary support.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "immersedBoundaryFvPatch.H"
#include "immersedBoundaryAdjustPhi.H"
#include "simpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::validOptions.insert("writep", "");
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
simpleControl simple(mesh);
# include "createIbMasks.H"
# include "createFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< nl << "Calculating potential flow" << endl;
// Do correctors over the complete set
while (simple.correctNonOrthogonal())
{
phi = faceIbMask*(linearInterpolate(U) & mesh.Sf());
// Adjust immersed boundary fluxes
immersedBoundaryAdjustPhi(phi, U);
// Adjust fluxes
adjustPhi(phi, U, p);
p.storePrevIter();
fvScalarMatrix pEqn
(
fvm::laplacian
(
dimensionedScalar
(
"1",
dimTime/p.dimensions()*dimensionSet(0, 2, -2, 0, 0),
1
),
p
)
==
fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
// Correct the flux
phi -= pEqn.flux();
if (!simple.finalNonOrthogonalIter())
{
p.relax();
}
Info<< "p min " << gMin(p.internalField())
<< " max " << gMax(p.internalField())
<< " masked min "
<< gMin(cellIbMask.internalField()*p.internalField())
<< " max "
<< gMax(cellIbMask.internalField()*p.internalField())
<< endl;
Info<< "continuity error = "
<< mag
(
fvc::div(faceIbMask*phi)
)().weightedAverage(mesh.V()).value()
<< endl;
Info<< "Contour continuity error = "
<< mag(sum(phi.boundaryField()))
<< endl;
U = fvc::reconstruct(phi);
U.correctBoundaryConditions();
Info<< "Interpolated U error = "
<< (
sqrt
(
sum
(
sqr
(
faceIbMask*
(
fvc::interpolate(U) & mesh.Sf()
)
- phi
)
)
)/sum(mesh.magSf())
).value()
<< endl;
}
// Calculate velocity magnitude
{
volScalarField magU = cellIbMask*mag(U);
Info << "IB-masked mag(U): max: " << gMax(magU.internalField())
<< " min: " << gMin(magU.internalField()) << endl;
}
// Force the write
U.write();
phi.write();
cellIbMask.write();
cellIbMaskExt.write();
if (args.optionFound("writep"))
{
// Find reference patch
label refPatch = -1;
scalar maxMagU = 0;
// Go through all velocity patches and find the one that fixes
// velocity to the largest value
forAll (U.boundaryField(), patchI)
{
const fvPatchVectorField& Upatch = U.boundaryField()[patchI];
if (Upatch.fixesValue())
{
// Calculate mean velocity
scalar u = sum(mag(Upatch));
label patchSize = Upatch.size();
reduce(u, sumOp<scalar>());
reduce(patchSize, sumOp<label>());
if (patchSize > 0)
{
scalar curMag = u/patchSize;
if (curMag > maxMagU)
{
refPatch = patchI;
maxMagU = curMag;
}
}
}
}
if (refPatch > -1)
{
// Calculate reference pressure
const fvPatchVectorField& Upatch = U.boundaryField()[refPatch];
const fvPatchScalarField& pPatch = p.boundaryField()[refPatch];
scalar patchE = sum(mag(pPatch + 0.5*magSqr(Upatch)));
label patchSize = Upatch.size();
reduce(patchE, sumOp<scalar>());
reduce(patchSize, sumOp<label>());
scalar e = patchE/patchSize;
Info<< "Using reference patch " << refPatch
<< " with mag(U) = " << maxMagU
<< " p + 0.5*U^2 = " << e << endl;
p.internalField() = e - 0.5*magSqr(U.internalField());
p.correctBoundaryConditions();
}
else
{
Info<< "No reference patch found. Writing potential function"
<< endl;
}
p.write();
}
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

3
applications/solvers/immersedBoundary/simpleIbFoam/Make/files
View file

@ -1,3 +0,0 @@
simpleIbFoam.C
EXE = $(FOAM_APPBIN)/simpleIbFoam

20
applications/solvers/immersedBoundary/simpleIbFoam/Make/options
View file

@ -1,20 +0,0 @@
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lfiniteVolume \
-lmeshTools \
-lsurfMesh \
-lsampling \
-ldynamicMesh \
-limmersedBoundary \
-limmersedBoundaryTurbulence \
-llduSolvers

10
applications/solvers/immersedBoundary/simpleIbFoam/UEqn.H
View file

@ -1,10 +0,0 @@
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff()
);
UEqn().relax();
solve(UEqn() == -fvc::grad(p));

42
applications/solvers/immersedBoundary/simpleIbFoam/createFields.H
View file

@ -1,42 +0,0 @@
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
mesh.schemesDict().setFluxRequired(p.name());
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, laminarTransport)
);

41
applications/solvers/immersedBoundary/simpleIbFoam/pEqn.H
View file

@ -1,41 +0,0 @@
{
volScalarField AU = UEqn().A();
U = UEqn().H()/AU;
// Immersed boundary update
U.correctBoundaryConditions();
UEqn.clear();
phi = faceIbMask*(fvc::interpolate(U) & mesh.Sf());
// Adjust immersed boundary fluxes
immersedBoundaryAdjustPhi(phi, U);
adjustPhi(phi, U, p);
// Non-orthogonal pressure corrector loop
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(1.0/AU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi -= pEqn.flux();
}
}
# include "immersedBoundaryContinuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
// Momentum corrector
U -= fvc::grad(p)/AU;
U.correctBoundaryConditions();
}

3
applications/solvers/incompressible/RichardsFoam/Make/files Executable file
View file

@ -0,0 +1,3 @@
RichardsFoam.C
EXE = $(FOAM_USER_APPBIN)/RichardsFoam

5
applications/solvers/incompressible/RichardsFoam/Make/options Executable file
View file

@ -0,0 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume

145
applications/solvers/incompressible/RichardsFoam/RichardsFoam.C Normal file
View file

@ -0,0 +1,145 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
RichardsFoam
Description
Transient solver for flow in unsaturated porous media
With chord slope formulation of the Richards equation.
van Genuchten laws for unsaturated hydraulic properties parametrisation
Global computation of the convergence criterium
Adaptative time stepping with a stabilisation procedure
NB 1: use backward scheme for time discretisation
NB 2: use only mesh with constant cell volumes
References
version 0.0 (develloped with OpenFOAM 2.0.1)
Details may be found in:
Orgogozo, L., Renon, N., Soulaine, C., Hénon, F., Tomer, S.K., Labat, D.,
Pokrovsky, O.S., Sekhar, M., Ababou, R., Quintard, M., Submitted.
Mechanistic modelling of water fluxes at the watershed scale: An open source
massively parallel solver for Richards equation.
Submitted to Computer Physics Communications.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "pimpleControl.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
// pimpleControl pimple(mesh);
# include "readPicardControls.H"
# include "createFields.H"
# include "initContinuityErrs.H"
# include "createTimeControls.H"
Info<< "\nStarting time loop\n" << endl;
// starting of the time loop.
while (runTime.loop())
{
// time step control operations.
# include "readTimeControls.H"
# include "setDeltaT.H"
// runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// Beginning of the stabilisation loop for the stabilised adaptive time
// step procedure.
for (int cyc = 0; cyc < nMaxCycle; cyc++)
{
// Beginning of the Picard loop.
for (int pic = 0; pic < nIterPicard; pic++)
{
# include "psiEqn.H"
}
// Exit test for the loop associated with the stabilisation cycles
// for the adaptive time step procedure.
if (crit < precPicard)
{
break;
}
else
{
Info << "Criterion not reached, restart time loop iteration"
<< "with a smaller time step / Error = " << crit
<< nl << endl;
runTime.setDeltaT((1/tFact)*runTime.deltaTValue());
Info<< "deltaT = " << runTime.deltaTValue() << endl;
}
// End of the stabilisation cycles loop.
}
// Warning test in case of convergence failure of the Picard loop.
if (crit >= precPicard)
{
Info<< "Convergence failure / Error = " << crit << nl << endl;
currentPicard = nIterPicard;
}
// Final updating of the result fields before going to the next time
// iteration.
psi_tmp = psi;
thtil_tmp = 0.5*
(
(1 + sign(psi_tmp)) + (1 - sign(psi_tmp))*
pow((1 + pow(mag(alpha*psi_tmp),n)), - (1 - (1/n)))
);
theta = (thetas - thetar)*thtil + thetar;
U = - Krel*((fvc::grad(psi)) + vuz);
// Writting of the result.
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
// end of the time loop.
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

209
applications/solvers/incompressible/RichardsFoam/createFields.H Normal file
View file

@ -0,0 +1,209 @@
Info<< "Reading transportProperties\n" << endl;
// reading of the physical constant associated with the considered problem.
// Localisation of the data within the considered case:
// constant/transportProperties
IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
);
dimensionedScalar K
(
transportProperties.lookup("K")
);
dimensionedScalar alpha
(
transportProperties.lookup("alpha")
);
dimensionedScalar thetas
(
transportProperties.lookup("thetas")
);
dimensionedScalar thetar
(
transportProperties.lookup("thetar")
);
dimensionedScalar n
(
transportProperties.lookup("n")
);
dimensionedScalar C
(
transportProperties.lookup("C")
);
dimensionedScalar S
(
transportProperties.lookup("S")
);
// declaration of the variable and results fields
// Water velocity field [m/s]
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Water saturation field [-]
Info<< "Reading field theta\n" << endl;
volScalarField theta
(
IOobject
(
"theta",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Water pressure field [m] - field of resolution.
Info<< "Reading field psi\n" << endl;
volScalarField psi
(
IOobject
(
"psi",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Field of residuals for the Picard loop [m].
Info<< "Reading field err\n" << endl;
volScalarField err
(
IOobject
(
"err",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Dimensionned unit scalar field [m].
Info<< "Reading field vuz\n" << endl;
volVectorField vuz
(
IOobject
(
"vuz",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Dimensionless unit vertical upward vector field.
Info<< "Reading field usf\n" << endl;
volScalarField usf
(
IOobject
(
"usf",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
# include "createPhi.H"
// Initialisation of the scalar containing the number of mesh cells. Note the
// use of gSum instead of sum.
double nbMesh;
nbMesh = gSum(usf);
// Initialisation of the scalar containing the residual for the exit test of the
// Picard loop.
double crit;
crit=0.;
// Initialisation of the token which counts the number of Picard iteraion for
// the adaptive time step procedure.
int currentPicard;
currentPicard = nIterPicard-3;
// Initialisation of the token which counts the number of Stabilisation cycle
// for the stabilisation of the adaptive time step procedure.
int sc;
sc = 0;
// Initialisation of the field of altitudes.
volVectorField positionVector = mesh.C();
volScalarField z = positionVector.component(vector::Z);
// Initialisation of the intermediate fields for the Picard loop.
volScalarField psi_tmp = psi;
volScalarField psim1 = psi;
// Initialisation of the varying transport properties for the Picard loop.
volScalarField thtil =
0.5*
(
(1 + sign(psi)) + (1 - sign(psi))*
pow((1 + pow(mag(alpha*psi),n)), - (1 - (1/n)))
);
volScalarField thtil_tmp =
0.5*
(
(1 + sign(psi_tmp)) + (1-sign(psi_tmp))*
pow((1 + pow(mag(alpha*psi_tmp),n)), - (1 - (1/n)))
);
volScalarField Krel =
0.5*
(
(1 + sign(psi))*K + (1 - sign(psi))*K*pow(thtil,0.5)*
pow((1 - pow((1 - pow(thtil,(n/(n - 1)))),(1 - (1/n)))),2)
);
volScalarField Crel =
S + 0.5*
(
(1 - sign(psi))*((thetas - thetar)*(thtil - thtil_tmp)*
(1./((usf*pos(psi - psi_tmp)*pos(psi_tmp - psi)) + psi - psi_tmp)))
);
// Initialisation of the gravity term.
volVectorField gradk = fvc::grad(Krel);
volScalarField gradkz = gradk.component(vector::Z);
// Initialisation of the velocity field.
U = - Krel*((fvc::grad(psi)) + vuz);

60
applications/solvers/incompressible/RichardsFoam/psiEqn.H Normal file
View file

@ -0,0 +1,60 @@
psim1 = psi;
psi = psi_tmp;
// Resolution of the linear system.
fvScalarMatrix psiEqn
(
Crel*fvm::ddt(psi)
==
fvm::laplacian(Krel, psi, "laplacian(Krel,psi)")
+ gradkz
);
psiEqn.relax();
psiEqn.solve();
// Update of the varying transport properties.
thtil = 0.5*
(
(1 + sign(psi)) + (1 - sign(psi))*
pow((1+pow(mag(alpha*psi),n)),-(1-(1/n)))
);
Krel = 0.5*
(
(1 + sign(psi))*K + (1 - sign(psi))*K*pow(thtil, 0.5)*
pow((1 - pow((1 - pow(thtil, (n/(n - 1)))),(1 - (1/n)))), 2)
);
Crel= S + 0.5*
(
(1 - sign(psi))*
(
(thetas - thetar)*(thtil - thtil_tmp)*
(
1./((usf*pos(psi - psi_tmp)*pos(psi_tmp - psi))
+ psi - psi_tmp)
)
)
);
// Update of the gravity term.
gradk = fvc::grad(Krel);
gradkz = gradk.component(2);
// Computation of the field of residuals for the exit test of
// the Picard loop.
err = psi - psim1;
// Computation of the residual for the exit test of the Picard
// loop. Note the use of gSum instead of sum.
crit = gSumMag(err)/nbMesh;
// exit test for the Picard loop.
if (crit < precPicard)
{
Info << " Erreur = " << crit
<< " Picard = " << pic
<< nl << endl;
currentPicard=pic;
break;
}

25
applications/solvers/incompressible/RichardsFoam/readPicardControls.H Normal file
View file

@ -0,0 +1,25 @@
// Reading of the controling parameters of the Picard loop and of the
// Stabilisation cycles loop for the adaptive time step procedure.
// Localisation of the data within the considered case: system/fvSolution
const dictionary& picardDict = mesh.solutionDict().subDict("Picard");
// Maximum number of Picard iterations.
const int nIterPicard =
picardDict.lookupOrDefault<int>("nIterPicard", 1);
// Maximum number of stabilisation cycles.
const int nMaxCycle =
picardDict.lookupOrDefault<int>("nMaxCycle", 1);
// Number of time iterations with low number of Picard iteration beyond which
// the time step is increased.
const int stabilisationThreshold =
picardDict.lookupOrDefault<int>("stabilisationThreshold", 1);
// Exit criterion for the Picard loop.
const double precPicard =
picardDict.lookupOrDefault<double>("precPicard", 1.);
// Factor of increase/decrease of the time step.
const double tFact =
picardDict.lookupOrDefault<double>("tFact", 1.);

82
applications/solvers/incompressible/RichardsFoam/setDeltaT.H Normal file
View file

@ -0,0 +1,82 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Global
setDeltaT
Description
Reset the timestep to maintain a low number of Picard iteration and a good
convergence. Reduction of time-step is immediate, but increase is damped to
avoid unstable oscillations.
\*---------------------------------------------------------------------------*/
scalar deltaTFact = 1.;
if (adjustTimeStep)
{
deltaTFact = 1.;
// If There is a too high number of Picard iteration, decrease of the time
// step.
if (currentPicard > nIterPicard-2)
{
deltaTFact = 1./tFact;
}
if (currentPicard >= 3)
{
sc = 0;
}
if ((currentPicard < 3) && (sc < stabilisationThreshold))
{
sc = sc+1;
}
// If there is more than 'stabilisationThreshold' iteration that the number
// of Picard iteration is low, increase of the time step.
if ((currentPicard < 3) && (sc == stabilisationThreshold))
{
deltaTFact = tFact;
sc = 0;
}
// Reset of the time step if needed (time step is always lower than
// 'maxDeltaT').
if (deltaTFact != 1.)
{
runTime.setDeltaT
(
min
(
deltaTFact*runTime.deltaTValue(),
maxDeltaT
)
);
}
Info<< "deltaT = " << runTime.deltaTValue() << endl;
}
// ************************************************************************* //

3
applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/files Executable file
View file

@ -0,0 +1,3 @@
spatialMeanValueRichardsonFoam.C
EXE = $(FOAM_USER_APPBIN)/spatialMeanValueRichardsonFoam

5
applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/options Executable file
View file

@ -0,0 +1,5 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume

110
applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/spatialMeanValueRichardsonFoam.C Normal file
View file

@ -0,0 +1,110 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
spatialMeanValue
Description
Calculates the spatial mean of the specified volScalarField over the
whole domain (for regular grid).
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "cyclicPolyPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
# include "addRegionOption.H"
timeSelector::addOptions();
argList::validArgs.append("fieldName");
# include "setRootCase.H"
# include "createTime.H"
instantList timeDirs = timeSelector::select0(runTime, args);
# include "createNamedMesh.H"
word fieldName(args.additionalArgs()[0]);
Info << "Calculating avarage pressure heads:" << endl;
forAll(timeDirs, timeI)
{
runTime.setTime(timeDirs[timeI], timeI);
IOobject fieldHeader
(
fieldName,
runTime.timeName(),
mesh,
IOobject::MUST_READ
);
// Check field exists
if (fieldHeader.headerOk())
{
mesh.readUpdate();
// Read field and calc mean, for regular grid
if (fieldHeader.headerClassName() == volScalarField::typeName)
{
volScalarField field(fieldHeader, mesh);
int nbMesh;
nbMesh = 0;
forAll(field, cellI)
{
nbMesh++;
}
Info<< runTime.timeName()<< " "
<< sum(field).value()/nbMesh<< " "
<< endl;
}
else
{
FatalError
<< "Only possible for volScalarField "<< "s "
<< nl << exit(FatalError);
}
}
else
{
Info<< " No field " << fieldName << endl;
}
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

32
applications/solvers/incompressible/pimpleDyMFoam/correctPhi.H
View file

@ -1,34 +1,6 @@
{
# include "continuityErrs.H"
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
for (label i=0; i<p.boundaryField().size(); i++)
{
if (p.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", p.dimensions(), 0),
pcorrTypes
);
volScalarField pcorr("pcorr", p);
pcorr *= 0;
// Initialise flux with interpolated velocity
phi = fvc::interpolate(U) & mesh.Sf();

4
applications/solvers/multiphase/interPhaseChangeFoam/createFields.H
View file

@ -104,7 +104,6 @@
label pdRefCell = 0;
scalar pdRefValue = 0.0;
setRefCell(pd, pimple.dict(), pdRefCell, pdRefValue);
mesh.schemesDict().setFluxRequired(pd.name());
scalar pRefValue = 0.0;
@ -128,3 +127,6 @@
(
incompressible::turbulenceModel::New(U, phi, twoPhaseProperties())
);
mesh.schemesDict().setFluxRequired(pd.name());
mesh.schemesDict().setFluxRequired(alpha1.name());

12
applications/utilities/immersedBoundary/refineImmersedBoundaryMesh/Make/options
View file

@ -1,13 +1,13 @@
EXE_INC = \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundaryDynamicMesh/lnInclude \
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lsurfMesh \
-lsampling \
-ldynamicMesh \
-limmersedBoundary
-limmersedBoundary \
-limmersedBoundaryDynamicMesh

7
applications/utilities/immersedBoundary/refineImmersedBoundaryMesh/refineImmersedBoundaryMesh.C
View file

@ -63,13 +63,6 @@ int main(int argc, char *argv[])
refineImmersedBoundaryMesh::IB_CELL_CELLS
);
}
else if (args.optionFound("ibCellCellFaces"))
{
rc = rib.refinementCells
(
refineImmersedBoundaryMesh::IB_CELL_CELL_FACES
);
}
else
{
FatalErrorIn(args.executable())

2
applications/utilities/immersedBoundary/writeIbMasks/Make/options
View file

@ -1,9 +1,11 @@
EXE_INC = \
-I$(LIB_SRC)/postProcessing/postCalc \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
EXE_LIBS = \
$(FOAM_LIBBIN)/postCalc.o \
-lfiniteVolume \
-lmeshTools \
-lsurfMesh \

126
applications/utilities/immersedBoundary/writeIbMasks/writeIbMasks.C
View file

@ -29,32 +29,128 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "calc.H"
#include "fvc.H"
#include "fvMatrices.H"
#include "immersedBoundaryFvPatch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createIbMasks.H"
Info<< nl << "Calculating gamma" << endl;
volScalarField gamma
(
IOobject
(
"gamma",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("one", dimless, 1)
);
gamma.internalField() = mesh.V()/mesh.cellVolumes();
gamma.write();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Report minimal live cell volume
scalar minLiveGamma = GREAT;
label minLiveCell = -1;
const scalarField& gammaIn = gamma.internalField();
Info<< "Time = " << runTime.timeName() << nl << endl;
forAll (mesh.boundary(), patchI)
{
if (isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
const immersedBoundaryFvPatch& ibPatch =
refCast<const immersedBoundaryFvPatch>
(
mesh.boundary()[patchI]
);
cellIbMask.write();
cellIbMaskExt.write();
const labelList& ibCells = ibPatch.ibPolyPatch().ibCells();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
forAll (ibCells, dcI)
{
if (gammaIn[ibCells[dcI]] < minLiveGamma)
{
minLiveGamma = gammaIn[ibCells[dcI]];
minLiveCell = ibCells[dcI];
}
}
}
}
Info<< "End\n" << endl;
Info<< "Min live cell " << minLiveCell
<< " gamma = " << minLiveGamma
<< endl;
return 0;
Info<< nl << "Calculating sGamma" << endl;
surfaceScalarField sGamma
(
IOobject
(
"sGamma",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("one", dimless, 0)
);
const surfaceScalarField& magSf = mesh.magSf();
const scalarField magFaceAreas = mag(mesh.faceAreas());
sGamma.internalField() =
magSf.internalField()/
scalarField::subField(magFaceAreas, mesh.nInternalFaces());
forAll (mesh.boundary(), patchI)
{
if (!isA<immersedBoundaryFvPatch>(mesh.boundary()[patchI]))
{
sGamma.boundaryField()[patchI] =
magSf.boundaryField()[patchI]/
mesh.boundary()[patchI].patchSlice(magFaceAreas);
gamma.boundaryField()[patchI] =
sGamma.boundaryField()[patchI];
}
}
sGamma.write();
// Check consistency of face area vectors
Info<< nl << "Calculating divSf" << endl;
volVectorField divSf
(
"divSf",
fvc::div(mesh.Sf())
);
divSf.write();
// Check divergence of face area vectors
scalarField magDivSf = mag(divSf)().internalField();
Info<< "Face areas divergence (min, max, average): "
<< "(" << min(magDivSf) << " " << max(magDivSf)
<< " " << average(magDivSf) << ")"
<< endl;
if (max(magDivSf) > 1e-9)
{
WarningIn("writeIbMasks")
<< "Possible problem with immersed boundary face area vectors: "
<< max(magDivSf)
<< endl;
}
Info<< endl;
}

3
applications/utilities/mesh/advanced/runDynamicMesh/Make/files Normal file
View file

@ -0,0 +1,3 @@
runDynamicMesh.C
EXE = $(FOAM_APPBIN)/runDynamicMesh

13
applications/solvers/immersedBoundary/icoDyMIbFoam/Make/options → applications/utilities/mesh/advanced/runDynamicMesh/Make/options
View file

@ -1,15 +1,12 @@
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
-I./lnInclude
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-limmersedBoundary \
-ldynamicFvMesh \
-ldynamicMesh \
-llduSolvers
-lmeshTools \
-lfiniteVolume \
-L$(MESQUITE_LIB_DIR) -lmesquite

13
applications/utilities/mesh/advanced/runDynamicMesh/createFields.H Normal file
View file

@ -0,0 +1,13 @@
Info<< "Reading field alpha\n" << endl;
volScalarField alpha
(
IOobject
(
"alpha",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

33
applications/solvers/immersedBoundary/simpleIbFoam/simpleIbFoam.C → applications/utilities/mesh/advanced/runDynamicMesh/runDynamicMesh.C
View file

@ -22,38 +22,33 @@ License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Application
simpleIbFoam
runDynamicMesh
Description
Steady-state solver for incompressible, turbulent flow
with immersed boundary support.
Creates dynamicFvMesh and updates it in a time loop. Reads indicator field
alpha used for dynamic mesh refinement to trigger refinement in certain
areas.
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved
Vuko Vukcevic, Wikki Ltd. All rights reserved.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
#include "dynamicFvMesh.H"
#include "simpleControl.H"
#include "immersedBoundaryFvPatch.H"
#include "immersedBoundaryAdjustPhi.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
# include "setRootCase.H"
# include "createTime.H"
# include "createMesh.H"
# include "createDynamicFvMesh.H"
simpleControl simple(mesh);
# include "createIbMasks.H"
# include "createFields.H"
# include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -63,15 +58,13 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << nl << endl;
// Pressure-velocity SIMPLE corrector
bool meshChanged = mesh.update();
// Write only if mesh has changed
if (meshChanged)
{
# include "UEqn.H"
# include "pEqn.H"
}
turbulence->correct();
runTime.write();
}
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
@ -80,7 +73,7 @@ int main(int argc, char *argv[])
Info<< "End\n" << endl;
return 0;
return(0);
}

1
applications/utilities/mesh/conversion/star4ToFoam/cellZoneToCellTableId/cellZoneToCellTableId.C
View file

@ -27,6 +27,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "objectRegistry.H"
#include "foamTime.H"
#include "volFields.H"
#include "OFstream.H"

1
applications/utilities/mesh/conversion/star4ToFoam/foamMeshToStar/foamMeshToStar.C
View file

@ -35,6 +35,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "objectRegistry.H"
#include "foamTime.H"
#include "volFields.H"
#include "cellModeller.H"

1
applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/files
View file

@ -4,6 +4,7 @@ linearNormal/linearNormal.C
linearRadial/linearRadial.C
sigmaRadial/sigmaRadial.C
wedge/wedge.C
gradedNormal/gradedNormal.C
LIB = $(FOAM_LIBBIN)/libextrudeModel

7
applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/options
View file

@ -1,8 +1,9 @@
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/ODE/lnInclude
EXE_LIBS = \
-lmeshTools \
-ldynamicMesh
-ldynamicMesh \
-lODE

143
applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.C Normal file
View file

@ -0,0 +1,143 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "gradedNormal.H"
#include "addToRunTimeSelectionTable.H"
#include "BisectionRoot.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace extrudeModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(gradedNormal, 0);
addToRunTimeSelectionTable(extrudeModel, gradedNormal, dictionary);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
gradedNormal::gradedNormal(const dictionary& dict)
:
extrudeModel(typeName, dict),
thickness_(readScalar(coeffDict_.lookup("thickness"))),
delta0_(readScalar(coeffDict_.lookup("initialCellLength"))),
expansionFactor_(1.0)
{
// Sanity checks
if (thickness_ <= SMALL)
{
FatalErrorIn("gradedNormal(const dictionary&)")
<< "thickness should be positive: " << thickness_
<< exit(FatalError);
}
if (delta0_ <= SMALL)
{
FatalErrorIn("gradedNormal(const dictionary&)")
<< "initialCellLength should be positive: " << delta0_
<< exit(FatalError);
}
const scalar maxExpFactor =
coeffDict_.lookupOrDefault<scalar>("maxExpansionFactor", 3.0);
if (maxExpFactor <= SMALL)
{
FatalErrorIn("gradedNormal(const dictionary&)")
<< "maxExpansionFactor should be positive: " << maxExpFactor
<< exit(FatalError);
}
const scalar bisectionTol =
coeffDict_.lookupOrDefault<scalar>("bisectionTol", 1e-5);
if (bisectionTol <= SMALL)
{
FatalErrorIn("gradedNormal(const dictionary&)")
<< "bisectionTolerance should be positive: " << bisectionTol
<< exit(FatalError);
}
// Create expansion factor equation represented as a function object
expansionFactorEqn eqn(*this);
// Calculate the expansionFactor using the bisection algorithm with the
// default tolerance of 1e-5
BisectionRoot<expansionFactorEqn> rootFinder
(
eqn,
bisectionTol
);
// Search for the root in [0, 3], where upper bound 3 is default
expansionFactor_ = rootFinder.root
(
0.0,
maxExpFactor
);
// Report the result
Info<< "Calculated expansion factor: " << expansionFactor_ << endl;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
gradedNormal::~gradedNormal()
{}
// * * * * * * * * * * * * * * * * Operators * * * * * * * * * * * * * * * * //
point gradedNormal::operator()
(
const point& surfacePoint,
const vector& surfaceNormal,
const label layer
) const
{
scalar d = 0.0;
for (label i = 0; i < layer; ++i)
{
d += delta0_*pow(expansionFactor_, i);
}
return surfacePoint + d*surfaceNormal;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace extrudeModels
} // End namespace Foam
// ************************************************************************* //

158
applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.H Normal file
View file

@ -0,0 +1,158 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::extrudeModels::gradedNormal
Description
Extrudes by transforming points normal to the surface. Input parameters are:
1. Extrusion thickness (in meters),
2. Initial delta (in meters),
3. Number of layers.
Expansion factor is calculated numerically using bisection algorithm.
Author
Vuko Vukcevic. FMENA Zagreb. All rights reserved.
\*---------------------------------------------------------------------------*/
#ifndef gradedNormal_H
#define gradedNormal_H
#include "point.H"
#include "extrudeModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace extrudeModels
{
/*---------------------------------------------------------------------------*\
Class gradedNormal Declaration
\*---------------------------------------------------------------------------*/
class gradedNormal
:
public extrudeModel
{
// Private data
//- Layer thickness
scalar thickness_;
//- Initial delta (cell size at the extruded patch)
scalar delta0_;
//- Expansion factor
scalar expansionFactor_;
// Expansion factor equation (functor class used by BisectionRoot)
class expansionFactorEqn
{
// Private data
//- Const reference to underlying gradedNormal model
const gradedNormal& gnm_;
public:
//- Construct given gradedNormal model
explicit expansionFactorEqn(const gradedNormal& gnm)
:
gnm_(gnm)
{}
//- Function call operator
scalar operator()(const scalar& x) const
{
scalar result = gnm_.thickness();
for (label i = 0; i <= gnm_.nLayers(); ++i)
{
result -= gnm_.delta0()*pow(x, i);
}
return result;
}
};
public:
//- Runtime type information
TypeName("gradedNormal");
// Constructors
//- Construct from components
gradedNormal(const dictionary& dict);
//- Destructor
~gradedNormal();
// Access functions
//- Return const reference to thickness
const scalar& thickness() const
{
return thickness_;
}
//- Return const reference to initial delta
const scalar& delta0() const
{
return delta0_;
}
// Member Operators
point operator()
(
const point& surfacePoint,
const vector& surfaceNormal,
const label layer
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace extrudeModels
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
applications/utilities/mesh/manipulation/setSet/setSet.C
View file

@ -894,7 +894,7 @@ int main(int argc, char *argv[])
Info<< "\nEnd\n" << endl;
return 0;
return !ok;
}

8
applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
View file

@ -660,7 +660,8 @@ autoPtr<mapPolyMesh> createRegionMesh
const EdgeMap<label>& interfaceToPatch,
const fvMesh& mesh,
const label regionI,
const word& regionName
const word& regionName,
const word& newMeshInstance
)
{
// Neighbour cellRegion.
@ -758,7 +759,7 @@ autoPtr<mapPolyMesh> createRegionMesh
IOobject
(
regionName,
mesh.time().timeName(),
newMeshInstance,
mesh.time(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
@ -977,7 +978,8 @@ void createAndWriteRegion
interfaceToPatch,
mesh,
regionI,
regionNames[regionI]
regionNames[regionI],
newMeshInstance
);
// Read in mesh as fvMesh for mapping. HJ, 19/May/2011

17
applications/utilities/parallelProcessing/decomposePar/decomposePar.C
View file

@ -284,6 +284,7 @@ int main(int argc, char *argv[])
(
"cellDist",
runTime.timeName(),
mesh.dbDir(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
@ -914,7 +915,7 @@ int main(int argc, char *argv[])
"faBoundary",
mesh.time().findInstance
(
mesh.dbDir()/polyMesh::meshSubDir,
mesh.meshDir(),
"boundary"
),
faMesh::meshSubDir,
@ -989,12 +990,20 @@ int main(int argc, char *argv[])
faMesh procMesh(procFvMesh);
word faceLabelsInstance =
procMesh.time().findInstance
(
procMesh.meshDir(),
"faceLabels"
);
labelIOList faceProcAddressing
(
IOobject
(
"faceProcAddressing",
"constant",
faceLabelsInstance,
procMesh.meshSubDir,
procFvMesh,
IOobject::MUST_READ,
@ -1007,7 +1016,7 @@ int main(int argc, char *argv[])
IOobject
(
"boundaryProcAddressing",
"constant",
faceLabelsInstance,
procMesh.meshSubDir,
procFvMesh,
IOobject::MUST_READ,
@ -1031,7 +1040,7 @@ int main(int argc, char *argv[])
IOobject
(
"edgeProcAddressing",
"constant",
faceLabelsInstance,
procMesh.meshSubDir,
procFvMesh,
IOobject::MUST_READ,

64
applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C
View file

@ -103,6 +103,37 @@ bool domainDecomposition::writeDecomposition()
label maxProcFaces = 0;
// Note: get cellLevel and pointLevel. Avoid checking whether they exist or
// not by hand. If they don't exist, simpy assume that the level is 0
const labelIOList globalCellLevel
(
IOobject
(
"cellLevel",
this->facesInstance(),
polyMesh::meshSubDir,
*this,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
labelList(nCells(), 0)
);
const labelIOList globalPointLevel
(
IOobject
(
"pointLevel",
this->facesInstance(),
polyMesh::meshSubDir,
*this,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
labelList(nPoints(), 0)
);
// Write out the meshes
for (label procI = 0; procI < nProcs_; procI++)
{
@ -615,6 +646,39 @@ bool domainDecomposition::writeDecomposition()
procBoundaryAddressing_[procI]
);
boundaryProcAddressing.write();
// Create and write cellLevel and pointLevel information
const unallocLabelList& cellMap = cellProcAddressing;
labelIOList procCellLevel
(
IOobject
(
"cellLevel",
procMesh.facesInstance(),
procMesh.meshSubDir,
procMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
labelList(globalCellLevel, cellMap)
);
procCellLevel.write();
const unallocLabelList& pointMap = pointProcAddressing;
labelIOList procPointLevel
(
IOobject
(
"pointLevel",
procMesh.facesInstance(),
procMesh.meshSubDir,
procMesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
labelList(globalPointLevel, pointMap)
);
procPointLevel.write();
}
Info<< nl

27
applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C
View file

@ -78,7 +78,7 @@ void faMeshDecomposition::distributeFaces()
IOobject
(
"faceProcAddressing",
"constant",
procMesh.facesInstance(),
procMesh.meshSubDir,
procMesh,
IOobject::MUST_READ,
@ -94,7 +94,7 @@ void faMeshDecomposition::distributeFaces()
IOobject
(
"owner",
"constant",
procMesh.facesInstance(),
procMesh.meshSubDir,
procMesh,
IOobject::MUST_READ,
@ -127,7 +127,7 @@ void faMeshDecomposition::distributeFaces()
IOobject
(
"faceProcAddressing",
"constant",
procMesh.facesInstance(),
procMesh.meshSubDir,
procMesh,
IOobject::MUST_READ,
@ -263,7 +263,7 @@ void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
IOobject
(
"pointProcAddressing",
"constant",
procFvMesh.facesInstance(),
procFvMesh.meshSubDir,
procFvMesh,
IOobject::MUST_READ,
@ -279,7 +279,7 @@ void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
IOobject
(
"faceProcAddressing",
"constant",
procFvMesh.facesInstance(),
procFvMesh.meshSubDir,
procFvMesh,
IOobject::MUST_READ,
@ -1202,7 +1202,7 @@ bool faMeshDecomposition::writeDecomposition()
IOobject
(
"boundaryProcAddressing",
"constant",
procFvMesh.facesInstance(),
procFvMesh.meshSubDir,
procFvMesh,
IOobject::MUST_READ,
@ -1340,13 +1340,20 @@ bool faMeshDecomposition::writeDecomposition()
maxProcPatches = max(maxProcPatches, nProcPatches);
maxProcEdges = max(maxProcEdges, nProcEdges);
word faceLabelsInstance =
procMesh.time().findInstance
(
procMesh.meshDir(),
"faceLabels"
);
// create and write the addressing information
labelIOList pointProcAddressing
(
IOobject
(
"pointProcAddressing",
"constant",
faceLabelsInstance,
procMesh.meshSubDir,
procFvMesh,
IOobject::NO_READ,
@ -1361,7 +1368,7 @@ bool faMeshDecomposition::writeDecomposition()
IOobject
(
"edgeProcAddressing",
"constant",
faceLabelsInstance,
procMesh.meshSubDir,
procFvMesh,
IOobject::NO_READ,
@ -1376,7 +1383,7 @@ bool faMeshDecomposition::writeDecomposition()
IOobject
(
"faceProcAddressing",
"constant",
faceLabelsInstance,
procMesh.meshSubDir,
procFvMesh,
IOobject::NO_READ,
@ -1391,7 +1398,7 @@ bool faMeshDecomposition::writeDecomposition()
IOobject
(
"boundaryProcAddressing",
"constant",
faceLabelsInstance,
procMesh.meshSubDir,
procFvMesh,
IOobject::NO_READ,

44
applications/utilities/parallelProcessing/decomposeSets/processorMeshes.C
View file

@ -157,6 +157,16 @@ Foam::polyMesh::readUpdateState Foam::processorMeshes::readUpdate()
// Check if any new meshes need to be read.
polyMesh::readUpdateState procStat = meshes_[procI].readUpdate();
// Force changed mesh if the processor mesh holds points. This is
// essential for changed region meshes and if the current time is
// equal to the initial time during construction.
// Pascal Beckstein, 14/Jul/2016
fileName timePath = meshes_[procI].time().timePath();
if (isFile(timePath/meshes_[procI].meshDir()/"points"))
{
procStat = polyMesh::POINTS_MOVED;
}
// Combine into overall mesh change status
if (stat == polyMesh::UNCHANGED)
{
@ -196,8 +206,36 @@ Foam::polyMesh::readUpdateState Foam::processorMeshes::readUpdate()
void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
{
// We may not use
// mesh.movePoints(newPoints);
// mesh.write();
// here, as this will fail for mesh regions containing
// directMappedPolyPatches or
// directMappedWallPolyPatches
// and probably other region-coupled patches. E.g. for both above
// mentioned directMapped patch types we have calls to
// directMappedPatchBase::map()
// in the corresponding initMovePoints(const pointField& p) member
// function. This however involves parallel communication and it would
// require all other connected regions to be already reconstructed and
// registered during reconstruction of the current region!
// Pascal Beckstein, 14/Jul/2016
// Create the new points
vectorField newPoints(mesh.nPoints());
pointIOField newPoints
(
IOobject
(
"points",
mesh.time().timeName(),
mesh.meshSubDir,
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh.allPoints()
);
forAll (meshes_, procI)
{
@ -225,8 +263,8 @@ void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
}
}
mesh.movePoints(newPoints);
mesh.write();
newPoints.write();
mesh.readUpdate();
}

44
applications/utilities/parallelProcessing/reconstructPar/processorMeshes.C
View file

@ -157,6 +157,16 @@ Foam::polyMesh::readUpdateState Foam::processorMeshes::readUpdate()
// Check if any new meshes need to be read.
polyMesh::readUpdateState procStat = meshes_[procI].readUpdate();
// Force changed mesh if the processor mesh holds points. This is
// essential for changed region meshes and if the current time is
// equal to the initial time during construction.
// Pascal Beckstein, 14/Jul/2016
fileName timePath = meshes_[procI].time().timePath();
if (isFile(timePath/meshes_[procI].meshDir()/"points"))
{
procStat = polyMesh::POINTS_MOVED;
}
// Combine into overall mesh change status
if (stat == polyMesh::UNCHANGED)
{
@ -196,8 +206,36 @@ Foam::polyMesh::readUpdateState Foam::processorMeshes::readUpdate()
void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
{
// We may not use
// mesh.movePoints(newPoints);
// mesh.write();
// here, as this will fail for mesh regions containing
// directMappedPolyPatches or
// directMappedWallPolyPatches
// and probably other region-coupled patches. E.g. for both above
// mentioned directMapped patch types we have calls to
// directMappedPatchBase::map()
// in the corresponding initMovePoints(const pointField& p) member
// function. This however involves parallel communication and it would
// require all other connected regions to be already reconstructed and
// registered during reconstruction of the current region!
// Pascal Beckstein, 14/Jul/2016
// Create the new points
vectorField newPoints(mesh.nPoints());
pointIOField newPoints
(
IOobject
(
"points",
mesh.time().timeName(),
mesh.meshSubDir,
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
mesh.allPoints()
);
forAll (meshes_, procI)
{
@ -225,8 +263,8 @@ void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
}
}
mesh.movePoints(newPoints);
mesh.write();
newPoints.write();
mesh.readUpdate();
}

14
applications/utilities/postProcessing/dataConversion/foamToTecplot360/Allwmake
View file

@ -5,10 +5,16 @@
cd ${0%/*} || exit 1 # run from this directory
set -x
# build tecio
wmake libso tecio/tecsrc
# foamToTecplot360 will not compile on MS Windows
if [ "$WM_OSTYPE" = "MSWindows" ]
then
echo MSWindows detected. Skipping compilation of foamToTecplot360.
else
# build tecio
wmake libso tecio/tecsrc
# build converter
wmake
# build converter
wmake
fi
# ----------------------------------------------------------------- end-of-file

2
applications/utilities/postProcessing/graphics/PV4FoamReader/Allwclean
View file

@ -3,6 +3,6 @@ cd ${0%/*} || exit 1 # run from this directory
set -x
rm -rf PV4FoamReader/Make
wclean libso vtkPV4Foam
wclean vtkPV4Foam
# ----------------------------------------------------------------- end-of-file

2
applications/utilities/postProcessing/wall/yPlusLES/yPlusLES.C
View file

@ -99,7 +99,7 @@ int main(int argc, char *argv[])
incompressible::LESModel::New(U, phi, laminarTransport)
);
volScalarField::GeometricBoundaryField d = nearWallDist(mesh).y();
volScalarField::GeometricBoundaryField d = nearWallDist(mesh);
volScalarField nuEff = sgsModel->nuEff();
const fvPatchList& patches = mesh.boundary();

2
applications/utilities/surface/surfaceCoarsen/Make/files
View file

@ -1,5 +1,5 @@
bunnylod/progmesh.C
bunnylod/vector.C
bunnylod/vectorb.C
surfaceCoarsen.C
EXE = $(FOAM_APPBIN)/surfaceCoarsen

2
applications/utilities/surface/surfaceCoarsen/bunnylod/bunnygut.C
View file

@ -19,7 +19,7 @@
#include <GL/gl.h>
#pragma warning(disable : 4244)
#include "vector.h"
#include "vectorb.h"
#include "font.h"
#include "progmesh.h"
#include "rabdata.h"

0
applications/utilities/surface/surfaceCoarsen/bunnylod/list.h → applications/utilities/surface/surfaceCoarsen/bunnylod/listb.h
View file

4
applications/utilities/surface/surfaceCoarsen/bunnylod/progmesh.C
View file

@ -13,8 +13,8 @@
#include <assert.h>
//#include <windows.h>
#include "vector.h"
#include "list.h"
#include "vectorb.h"
#include "listb.h"
#include "progmesh.h"
#define min(x,y) (((x) <= (y)) ? (x) : (y))

4
applications/utilities/surface/surfaceCoarsen/bunnylod/progmesh.h
View file

@ -18,8 +18,8 @@
#ifndef PROGRESSIVE_MESH_H
#define PROGRESSIVE_MESH_H
#include "vector.h"
#include "list.h"
#include "vectorb.h"
#include "listb.h"
class tridata {
public:

2
applications/utilities/surface/surfaceCoarsen/bunnylod/vector.C → applications/utilities/surface/surfaceCoarsen/bunnylod/vectorb.C
View file

@ -3,7 +3,7 @@
#include <math.h>
#include <assert.h>
#include "vector.h"
#include "vectorb.h"
float sqr(float a) {return a*a;}

0
applications/utilities/surface/surfaceCoarsen/bunnylod/vector.h → applications/utilities/surface/surfaceCoarsen/bunnylod/vectorb.h
View file

30
bin/tools/RunFunctions
View file

@ -73,6 +73,11 @@ runParallel ()
LOG_NAME=log.$APP_NAME
fi
if [ "$WM_OSTYPE" = "MSWindows" ]
then
APP_RUN="${APP_RUN}.exe"
fi
if [ -f $LOG_NAME ] ; then
echo "$APP_NAME already run on $PWD: remove log file to run"
else
@ -139,6 +144,17 @@ cloneCase ()
makeFsiCaseLinks ()
{
if [ "$WM_OSTYPE" = "MSWindows" ]
then
cd $1
cd system
cp -r ../../$2/system $2
cd ../constant
cp -r ../../$2/constant $2
cd ../0
cp -r ../../$2/0 $2
cd ../..
else
cd $1
cd system
ln -s ../../$2/system $2
@ -147,10 +163,23 @@ makeFsiCaseLinks ()
cd ../0
ln -s ../../$2/0 $2
cd ../..
fi
}
makeFsiResultsLinks ()
{
if [ "$WM_OSTYPE" = "MSWindows" ]
then
cd $1
TIME_DIRS=`foamInfoExec -times | sed '1,/constant/d'`
echo "makeFsiResultsLinks for" $TIME_DIRS
cd ../$2
for T in $TIME_DIRS
do
cp -r ../$1/${T}/solid ${T}
done
cd ..
else
cd $1
TIME_DIRS=`foamInfoExec -times | sed '1,/constant/d'`
echo "makeFsiResultsLinks for" $TIME_DIRS
@ -160,6 +189,7 @@ makeFsiResultsLinks ()
ln -s ../$1/${T}/solid ${T}
done
cd ..
fi
}
#------------------------------------------------------------------------------

24
etc/bashrc
View file

@ -160,10 +160,11 @@ done
: ${WM_OSTYPE:=POSIX}; export WM_OSTYPE
# Compiler: set to Gcc or Icc (for Intel's icc)
# Compiler: set to Gcc, Icc (for Intel's icc) or Clang
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
: ${WM_COMPILER:=Gcc}; export WM_COMPILER
#: ${WM_COMPILER:=Icc}; export WM_COMPILER
#: ${WM_COMPILER:=Clang}; export WM_COMPILER
export WM_COMPILER_ARCH=
export WM_COMPILER_LIB_ARCH=
@ -316,6 +317,19 @@ Linux)
export WM_LDFLAGS='-mfloat-abi=hard'
;;
aarch64)
# https://gcc.gnu.org/onlinedocs/gcc/AArch64-Options.html
# This command will provide more details for your specific processor
# gcc -c -Q -march=native -mtune=native --help=target
WM_ARCH=linuxARM8
export WM_COMPILER_LIB_ARCH=64
export WM_CC='gcc'
export WM_CXX='g++'
export WM_CFLAGS='-fPIC -march=native -mtune=native'
export WM_CXXFLAGS='-fPIC -march=native -mtune=native'
export WM_LDFLAGS=''
;;
*)
echo Unknown processor type `uname -m` for Linux
;;
@ -366,7 +380,7 @@ Darwin)
echo "Using Macports binaries"
fi
export WM_USE_MACPORT=1
export WM_USE_MACPORT=0
export WM_BASE_COMPILER=`echo $WM_COMPILER | tr -d "[:digit:]"`
export WM_MACPORT_MPI_VERSION=`echo $WM_COMPILER | tr "[:upper:]" "[:lower:]"`
export WM_MACPORT_VERSION=`echo $WM_MACPORT_MPI_VERSION | tr -d "[:alpha:]" | sed -e "s/\(.\)\(.\)/\1\.\2/"`
@ -434,9 +448,15 @@ Darwin)
export WM_CXX="g++-mp-$WM_MACPORT_VERSION"
export WM_FC="gfortran-mp-$WM_MACPORT_VERSION"
elif [ "$WM_BASE_COMPILER" == "Clang" ]
then
if [ "$WM_USE_MACPORT" == "1" ]
then
export WM_CC="clang-mp-$WM_MACPORT_VERSION"
export WM_CXX="clang++-mp-$WM_MACPORT_VERSION"
else
export WM_CC="clang"
export WM_CXX="clang++"
fi
# Seems like there is no Fortran-frontend for LLVM at thistime
elif [ "$WM_BASE_COMPILER" == "Dragonegg" ]
then

4
etc/prefs.sh-EXAMPLE
View file

@ -160,8 +160,10 @@ export FOAM_VERBOSE=1
# System installed bison
#export BISON_SYSTEM=1
# System installed flex
# System installed flex. FLEX_DIR should point to the directory where
# $FLEX_DIR/bin/flex is located
#export FLEX_SYSTEM=1
#export FLEX_DIR=/usr
# System installed m4
#export M4_SYSTEM=1

3
src/Allwmake
View file

@ -77,8 +77,7 @@ wmake libso solidModels
wmake libso dbns
wmake libso immersedBoundary/immersedBoundary
wmake libso immersedBoundary/immersedBoundaryTurbulence
wmake libso immersedBoundary/immersedBoundaryForce
#wmake libso immersedBoundary/immersedBoundaryTurbulence
wmake libso immersedBoundary/immersedBoundaryDynamicMesh
( cd cudaSolvers ; ./Allwmake )

19
src/ODE/Make/files
View file

@ -15,6 +15,23 @@ $(translationODE)/translationODE.C
sixDOF = sixDOF
$(sixDOF)/finiteRotation/finiteRotation.C
$(sixDOF)/sixDOFqODE/sixDOFqODE.C
$(sixDOF)/sixDOFbodies/sixDOFbodies.C
$(sixDOF)/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.C
$(sixDOF)/sixDOFODE/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.C
$(sixDOF)/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.C
$(sixDOF)/sixDOFODE/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.C
$(sixDOF)/sixDOFODE/constraints/translationalConstraints/periodicOscillation/periodicOscillation.C
$(sixDOF)/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.C
$(sixDOF)/sixDOFODE/restraints/translationalRestraints/linearSpringDamper/linearSpringDamper.C
$(sixDOF)/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.C
$(sixDOF)/sixDOFODE/restraints/rotationalRestraints/angularDamper/angularDamper.C
$(sixDOF)/sixDOFODE/sixDOFODE.C
$(sixDOF)/sixDOFODE/newSixDOFODE.C
$(sixDOF)/quaternionSixDOF/quaternionSixDOF.C
$(sixDOF)/geometricSixDOF/geometricSixDOF.C
$(sixDOF)/sixDOFBodies/sixDOFBodies.C
LIB = $(FOAM_LIBBIN)/libODE

6
src/ODE/ODESolvers/ODESolver/ODESolver.H
View file

@ -53,7 +53,7 @@ class ODESolver
protected:
// Private data
// Protected data
//- Reference to ODE
ODE& ode_;
@ -62,6 +62,8 @@ protected:
mutable scalarField dydx_;
private:
// Private Member Functions
//- Disallow default bitwise copy construct
@ -126,7 +128,7 @@ public:
) const = 0;
virtual void solve
void solve
(
const scalar xStart,
const scalar xEnd,

54
src/ODE/sixDOF/finiteRotation/finiteRotation.C
View file

@ -26,22 +26,21 @@ Class
Author
Dubravko Matijasevic, FSB Zagreb. All rights reserved.
Update by Hrvoje Jasak
Hrvoje Jasak, FSB Zagreb. All rights reserved.
Vuko Vukcevic, FSB Zagreb. All rights reserved.
\*---------------------------------------------------------------------------*/
#include "objectRegistry.H"
#include "finiteRotation.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Static Functions * * * * * * * * * * * * * //
Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
{
vector ur = - *( inv(I + rotT) & (I - rotT) );
// Scaling to a unit vector. HJ, problems with round-off
// HJ, 4/Aug/2008
// Scaling to a unit vector. Problems with round-off. HJ, 4/Aug/2008
if (mag(ur) > SMALL)
{
return ur/(mag(ur) + SMALL);
@ -49,8 +48,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
else
{
// Rotation vector is undertermined at zero rotation
// Returning arbitrary unit vector
// HJ, 4/Mar/2015
// Returning arbitrary unit vector. HJ, 4/Mar/2015
return vector(0, 0, 1);
}
}
@ -79,22 +77,15 @@ Foam::vector Foam::finiteRotation::eulerAngles(const tensor& rotT)
scalar& pitchAngle = eulerAngles.y();
scalar& yawAngle = eulerAngles.z();
// Calculate roll angle
rollAngle = atan2(rotT.yz(), rotT.zz());
// Use mag to avoid negative value due to round-off
// HJ, 24/Feb/2016
// Bugfix: sqr. SS, 18/Apr/2016
const scalar c2 = sqrt(sqr(rotT.xx()) + sqr(rotT.xy()));
// Calculate pitch angle
pitchAngle = atan2(-rotT.xz(), c2);
pitchAngle = asin(rotT.xz());
const scalar s1 = sin(rollAngle);
const scalar c1 = cos(rollAngle);
// Calculate roll angle
const scalar cosPitch = cos(pitchAngle);
rollAngle = asin(-rotT.yz()/cosPitch);
// Calculate yaw angle
yawAngle = atan2(s1*rotT.zx() - c1*rotT.yx(), c1*rotT.yy() - s1*rotT.zy());
yawAngle = asin(-rotT.xy()/cosPitch);
return eulerAngles;
}
@ -122,6 +113,14 @@ Foam::finiteRotation::finiteRotation
{}
Foam::finiteRotation::finiteRotation(const tensor& R)
:
eInitial_(R),
rotTensor_(R),
rotIncrementTensor_(tensor::zero)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::finiteRotation::~finiteRotation()
@ -130,12 +129,23 @@ Foam::finiteRotation::~finiteRotation()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::finiteRotation::updateRotation(const tensor& R)
{
rotIncrementTensor_ = (R & rotTensor_.T());
rotTensor_ = R;
}
void Foam::finiteRotation::updateRotation(const HamiltonRodriguezRot& rot)
{
tensor rotR = rot.R();
updateRotation(rot.R());
}
rotIncrementTensor_ = (rotR & (rotTensor_.T()));
rotTensor_ = rotR;
void Foam::finiteRotation::updateRotationWithIncrement(const tensor& incR)
{
rotIncrementTensor_ = incR;
rotTensor_ = (incR & rotTensor_);
}

42
src/ODE/sixDOF/finiteRotation/finiteRotation.H
View file

@ -26,7 +26,8 @@ Class
Author
Dubravko Matijasevic, FSB Zagreb. All rights reserved.
Updated by Hrvoje Jasak
Hrvoje Jasak, FSB Zagreb. All rights reserved.
Vuko Vukcevic, FSB Zagreb. All rights reserved.
SourceFiles
finiteRotation.C
@ -61,20 +62,6 @@ class finiteRotation
tensor rotIncrementTensor_;
// Private Member Functions
//- Calculate unit rotation vector from given rotation tensor
static vector rotVector(const tensor& rotT);
//- Calculate rotation angle from given rotation tensor
static scalar rotAngle(const tensor& rotT);
//- Calculate Euler angles (x-y-z (roll-pitch-yaw) convention) from
// given rotation tensor. Reference: Mike Day, Insomniac Games,
// Extracting Euler Angles from a Rotation Matrix.
static vector eulerAngles(const tensor& rotT);
public:
// Constructors
@ -89,17 +76,40 @@ public:
const scalar& angle
);
//- Construct from rotation tensor
explicit finiteRotation(const tensor& R);
// Destructor
~finiteRotation();
// Static Functions
//- Calculate unit rotation vector from given rotation tensor
static vector rotVector(const tensor& rotT);
//- Calculate rotation angle from given rotation tensor
static scalar rotAngle(const tensor& rotT);
//- Calculate Euler angles (x-y-z (roll-pitch-yaw) convention by Bryan)
// given the rotation tensor (global-to-local transformation).
// Reference: Nikravesh: Computer-Aided Analysis of Mechanical Systems
static vector eulerAngles(const tensor& rotT);
// Member Functions
//- Update rotation
//- Update rotation given rotation tensor
void updateRotation(const tensor& R);
//- Update rotation given HamiltonRodriguezRot (quaternions)
void updateRotation(const HamiltonRodriguezRot& rot);
//- Update rotation given increment rotation tensor
void updateRotationWithIncrement(const tensor& incR);
//- Return initial quaternions
const HamiltonRodriguezRot& eInitial() const;

646
src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C Normal file
View file

@ -0,0 +1,646 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
geometricSixDOF
\*---------------------------------------------------------------------------*/
#include "geometricSixDOF.H"
#include "scalarSquareMatrix.H"
#include "OutputControlDictionary.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(geometricSixDOF, 0);
addToRunTimeSelectionTable(sixDOFODE, geometricSixDOF, dictionary);
}
const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncTensorTol_
(
"geometricSixDOFRotIncTensorTol",
1e-10
);
const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncRateTol_
(
"geometricSixDOFRotIncRateTol",
1e-6
);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::dimensionedVector Foam::geometricSixDOF::A
(
const dimensionedVector& xR,
const dimensionedVector& uR,
const tensor& R,
const scalar t
) const
{
// Create a scalar square matrix representing Newton equations with
// constraints and the corresponding source (right hand side vector).
// Note: size of the matrix is 3 + number of constraints
scalarField rhs(translationalConstraints().size() + 3, 0.0);
scalarSquareMatrix M(rhs.size(), 0.0);
// Insert mass and explicit forcing into the system. Note: translations are
// solved in the global coordinate system and the explicit forcing contains
// restraining forces
const dimensionedVector explicitForcing =
force
(
t,
R.T(),
xR.value(),
uR.value()
);
const vector& efVal = explicitForcing.value();
const scalar& m = mass().value();
forAll(efVal, dirI)
{
M[dirI][dirI] = m;
rhs[dirI] = efVal[dirI];
}
// Insert contributions from the constraints
forAll(translationalConstraints(), tcI)
{
// Get reference to current constraint
const translationalConstraint& curTc = translationalConstraints()[tcI];
// Get matrix contribution from constraint
const vector mc =
curTc.matrixContribution
(
t,
R.T(),
xR.value(),
uR.value()
);
// Get matrix index
const label index = tcI + 3;
// Insert contributions into the matrix
forAll(mc, dirI)
{
M[dirI][index] = mc[dirI];
M[index][dirI] = mc[dirI];
}
// Insert source contribution (remainder of the constraint function)
rhs[index] = curTc.sourceContribution(t, R.T(), xR.value(), uR.value());
}
// Solve the matrix using LU decomposition. Note: solution is in the rhs and
// it contains accelerations in the first three entries and corresponding
// Lagrangian multipliers in other entries.
scalarSquareMatrix::LUsolve(M, rhs);
return
dimensionedVector
(
"A",
force().dimensions()/mass().dimensions(),
vector(rhs[0], rhs[1], rhs[2])
);
}
Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
(
const tensor& R,
const dimensionedVector& omega,
const scalar t
) const
{
// Create a scalar square matrix representing Euler equations with
// constraints and the corresponding source (right hand side vector).
// Note: size of the matrix is 3 + number of constraints
scalarField rhs(rotationalConstraints().size() + 3, 0.0);
scalarSquareMatrix J(rhs.size(), 0.0);
// Get current inertial-to-local transformation
const dimensionedTensor RT("RT", dimless, R.T());
// Insert moment of inertia and explicit forcing into the system
const dimensionedVector explicitForcing
(
E(omega) // Euler part
+ (
RT
& moment
(
t,
RT.value(),
omega.value()
)
) // External torque with restraints
);
const vector& efVal = explicitForcing.value();
const diagTensor& I = momentOfInertia().value();
forAll(efVal, dirI)
{
J[dirI][dirI] = I[dirI];
rhs[dirI] = efVal[dirI];
}
// Insert contributions from the constraints
forAll(rotationalConstraints(), rcI)
{
// Get reference to current constraint
const rotationalConstraint& curRc = rotationalConstraints()[rcI];
// Get matrix contribution from the constraint
const vector mc =
curRc.matrixContribution
(
t,
RT.value(),
omega.value()
);
// Get matrix index
const label index = rcI + 3;
// Insert contributions into the matrix
forAll(mc, dirI)
{
J[dirI][index] = mc[dirI];
J[index][dirI] = mc[dirI];
}
// Insert source contribution (remainder of the constraint function)
rhs[index] = curRc.sourceContribution(t, RT.value(), omega.value());
}
// Solve the matrix using LU decomposition. Note: solution is in the rhs and
// it contains OmegaDot's in the first three entries and corresponding
// Lagrangian multipliers in other entries.
scalarSquareMatrix::LUsolve(J, rhs);
return
dimensionedVector
(
"OmegaDot",
moment().dimensions()/momentOfInertia().dimensions(),
vector(rhs[0], rhs[1], rhs[2])
);
}
Foam::dimensionedVector Foam::geometricSixDOF::E
(
const dimensionedVector& omega
) const
{
return (*(momentOfInertia() & omega) & omega);
}
Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const
{
tensor R;
// Calculate the magnitude of the rotational increment vector
const scalar magRotInc = mag(rotInc);
if (magRotInc < rotIncTensorTol_())
{
// No rotational increment
R = I;
}
else
{
// Calculate rotational increment tensor using the exponential map
// Skew-symmetric tensor corresponding to the rotation increment
const tensor skewRotInc(*rotInc);
R = I
+ skewRotInc*sin(magRotInc)/magRotInc
+ (skewRotInc & skewRotInc)*(1.0 - cos(magRotInc))/sqr(magRotInc);
}
return R;
}
Foam::vector Foam::geometricSixDOF::dexpMap
(
const vector& rotInc,
const vector& omega
) const
{
vector rotIncDot;
// Calculate the magnitude of the rotational increment vector
const scalar magRotInc = mag(rotInc);
if (magRotInc < rotIncRateTol_())
{
// Stabilised calculation of rotation increment to avoid small
// denominators
const tensor lbRotIncOmega(lieBracket(rotInc, omega));
rotIncDot =
(
I
+ lbRotIncOmega/2.0
+ lieBracket(rotInc, *lbRotIncOmega)/12.0
) & omega;
}
else
{
// Calculate rate of the rotational increment vector using the
// differential of the exponential map
// Skew-symmetric tensor corresponding to the rotation increment
const tensor skewRotInc(*rotInc);
rotIncDot =
(
I
+ 0.5*skewRotInc
- (skewRotInc & skewRotInc)*
(magRotInc/tan(magRotInc/2.0) - 2.0)/(2.0*sqr(magRotInc))
) & omega;
}
return rotIncDot;
}
// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
{
// First set the state in base class subobject
sixDOFODE::setState(sd);
// Cast sixDOFODE& to geometricSixDOF&
const geometricSixDOF& gsd = refCast<const geometricSixDOF>(sd);
// Set state variables for this class
Xrel_ = gsd.Xrel_;
U_ = gsd.U_;
Uaverage_ = gsd.Uaverage_;
rotation_ = gsd.rotation_;
omega_ = gsd.omega_;
omegaAverage_ = gsd.omegaAverage_;
// Copy ODE coefficients: this carries actual ODE state
// HJ, 23/Mar/2015
coeffs_ = gsd.coeffs_;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
:
sixDOFODE(io),
Xrel_(dict().lookup("Xrel")),
U_(dict().lookup("U")),
Uaverage_("Uaverage", U_),
rotation_(tensor(dict().lookup("rotationTensor"))),
rotIncrement_
(
dict().lookupOrDefault<tensor>("rotationIncrementTensor", tensor::zero)
),
omega_(dict().lookup("omega")),
omegaAverage_("omegaAverage", omega_),
coeffs_(12, 0.0)
{
// Set ODE coefficients from position and rotation
// Linear displacement relative to spring equilibrium
const vector& Xval = Xrel_.value();
coeffs_[0] = Xval.x();
coeffs_[1] = Xval.y();
coeffs_[2] = Xval.z();
// Linear velocity
const vector& Uval = U_.value();
coeffs_[3] = Uval.x();
coeffs_[4] = Uval.y();
coeffs_[5] = Uval.z();
// Rotational velocity in non - inertial coordinate system
const vector& omegaVal = omega_.value();
coeffs_[6] = omegaVal.x();
coeffs_[7] = omegaVal.y();
coeffs_[8] = omegaVal.z();
// Increment of the rotation vector (zero for initial condition)
coeffs_[9] = 0;
coeffs_[10] = 0;
coeffs_[11] = 0;
}
Foam::geometricSixDOF::geometricSixDOF
(
const word& name,
const geometricSixDOF& gsd
)
:
sixDOFODE(name, gsd),
Xrel_(gsd.Xrel_.name(), gsd.Xrel_),
U_(gsd.U_.name(), gsd.U_),
Uaverage_(gsd.Uaverage_.name(), gsd.Uaverage_),
rotation_(gsd.rotation_),
omega_(gsd.omega_.name(), gsd.omega_),
omegaAverage_(gsd.omegaAverage_.name(), gsd.omegaAverage_),
coeffs_(gsd.coeffs_)
{}
Foam::autoPtr<Foam::sixDOFODE> Foam::geometricSixDOF::clone
(
const word& name
) const
{
// Create and return an autoPtr to the new geometricSixDOF object with a
// new name
return autoPtr<sixDOFODE>
(
new geometricSixDOF
(
name,
*this
)
);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::geometricSixDOF::~geometricSixDOF()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::dimensionedVector& Foam::geometricSixDOF::Xrel() const
{
return Xrel_;
}
const Foam::dimensionedVector& Foam::geometricSixDOF::omega() const
{
return omega_;
}
Foam::dimensionedVector Foam::geometricSixDOF::X() const
{
return Xequilibrium() + Xrel_;
}
const Foam::dimensionedVector& Foam::geometricSixDOF::U() const
{
return U_;
}
const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const
{
return Uaverage_;
}
const Foam::dimensionedVector& Foam::geometricSixDOF::omegaAverage() const
{
return omegaAverage_;
}
Foam::tensor Foam::geometricSixDOF::toRelative() const
{
return rotation_.T();
}
Foam::tensor Foam::geometricSixDOF::toAbsolute() const
{
return rotation_;
}
const Foam::tensor& Foam::geometricSixDOF::rotIncrementTensor() const
{
return rotIncrement_;
}
void Foam::geometricSixDOF::derivatives
(
const scalar x,
const scalarField& y,
scalarField& dydx
) const
{
// Translation
// Set the derivatives for displacement
dydx[0] = y[3];
dydx[1] = y[4];
dydx[2] = y[5];
dimensionedVector curX("curX", dimLength, vector(y[0], y[1], y[2]));
dimensionedVector curU("curU", dimVelocity, vector(y[3], y[4], y[5]));
// Get rotational increment vector (u)
const vector rotIncrementVector(y[9], y[10], y[11]);
// Calculate current rotation tensor obtained with exponential map
const tensor curRot = (rotation_ & expMap(rotIncrementVector));
// Calculate translational acceleration using current rotation
const vector accel = A(curX, curU, curRot, x).value();
// Set the derivatives for velocity
dydx[3] = accel.x();
dydx[4] = accel.y();
dydx[5] = accel.z();
// Rotation
dimensionedVector curOmega
(
"curOmega",
dimless/dimTime,
vector(y[6], y[7], y[8])
);
// Calculate rotational acceleration using current rotation
const vector omegaDot = OmegaDot(curRot, curOmega, x).value();
dydx[6] = omegaDot.x();
dydx[7] = omegaDot.y();
dydx[8] = omegaDot.z();
// Calculate derivative of rotIncrementVector using current rotation
// information
const vector rotIncrementVectorDot =
dexpMap(rotIncrementVector, curOmega.value());
// Set the derivatives for rotation
dydx[9] = rotIncrementVectorDot.x();
dydx[10] = rotIncrementVectorDot.y();
dydx[11] = rotIncrementVectorDot.z();
}
void Foam::geometricSixDOF::update(const scalar delta)
{
// Translation
// Get displacement
const vector Xold = Xrel_.value();
vector& Xval = Xrel_.value();
Xval.x() = coeffs_[0];
Xval.y() = coeffs_[1];
Xval.z() = coeffs_[2];
// Get velocity
vector& Uval = U_.value();
Uval.x() = coeffs_[3];
Uval.y() = coeffs_[4];
Uval.z() = coeffs_[5];
// Stabilise translational constraints if necessary
forAll(translationalConstraints(), tcI)
{
// Note: get end time for this ODE step from mesh, assuming that the top
// level solves this ode from [t - deltaT, t], yielding solution at
// t. Done this way to preserve the interface of ODE class.
// VV, 10/Mar/2017.
const scalar t = dict().time().value();
translationalConstraints()[tcI].stabilise(t, Xval, Uval);
}
// Update (possibly constrained) displacement
coeffs_[0] = Xval.x();
coeffs_[1] = Xval.y();
coeffs_[2] = Xval.z();
// Update (possibly constrained) velocity
coeffs_[3] = Uval.x();
coeffs_[4] = Uval.y();
coeffs_[5] = Uval.z();
// Update average velocity
Uaverage_.value() = (Xval - Xold)/delta;
// Rotation
// Get angular velocity
const vector omegaOld = omega_.value();
vector& omegaVal = omega_.value();
omegaVal.x() = coeffs_[6];
omegaVal.y() = coeffs_[7];
omegaVal.z() = coeffs_[8];
// Update rotational increment tensor
rotIncrement_ = expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]));
// Update rotational tensor
rotation_ = (rotation_ & rotIncrement_);
// Stabilise rotational constraints if necessary
forAll(rotationalConstraints(), rcI)
{
// Note: get end time for this ODE step from mesh, assuming that the top
// level solves this ode from [t - deltaT, t], yielding solution at
// t. Done this way to preserve the interface of ODE class.
// VV, 10/Mar/2017.
const scalar t = dict().time().value();
rotationalConstraints()[rcI].stabilise(t, omegaVal);
}
// Update (possibly constrained) omega
coeffs_[6] = omegaVal.x();
coeffs_[7] = omegaVal.y();
coeffs_[8] = omegaVal.z();
// Reset increment vector in coefficients for the next step
coeffs_[9] = 0;
coeffs_[10] = 0;
coeffs_[11] = 0;
// Update average omega
omegaAverage_.value() = 0.5*(omegaVal + omegaOld);
}
bool Foam::geometricSixDOF::writeData(Ostream& os) const
{
// First write the part related to base class subobject
sixDOFODE::writeData(os);
// Write type name
os.writeKeyword("type") << tab << type()
<< token::END_STATEMENT << nl << nl;
// Write data
os.writeKeyword("Xrel") << tab << Xrel_
<< token::END_STATEMENT << nl;
os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
os.writeKeyword("rotationTensor") << tab << rotation_
<< token::END_STATEMENT << nl;
os.writeKeyword("rotationIncrementTensor") << tab << rotIncrement_
<< token::END_STATEMENT << nl;
os.writeKeyword("omega") << tab << omega_
<< token::END_STATEMENT << nl << endl;
return os.good();
}
// ************************************************************************* //

319
src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H Normal file
View file

@ -0,0 +1,319 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
General Public License for more details.
You should have received a copy of the GNU General Public License
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
geometricSixDOF
Description
6-DOF solver using a geometric method for integration of rotations.
Run-time selectable constraints are handled via Lagrangian multipliers using
the interface from translational/rotationalConstraint classes.
Reference (bibtex):
@article {terzeEtAl2016,
Author = {Terze, Z. and M\"{u}ller, A. and Zlatar, D.},
title = {Singularity-free time integration of rotational quaternions
using non-redundant ordinary differential equations},
Journal = {Multibody System Dynamics},
Year = {2016},
Volume = {38},
Number = {3},
Pages = {201--225}
}
Note on convention: rotation tensor (R, or rotation_) defines
body-to-inertial coordinate system transformation
(local-to-global). Opposite as in quaternionSixDOF.
Author
Viktor Pandza, FSB Zagreb. All rights reserved.
Vuko Vukcevic, FSB Zagreb. All rights reserved.
SourceFiles
geometricSixDOF.C
\*---------------------------------------------------------------------------*/
#ifndef geometricSixDOF_H
#define geometricSixDOF_H
#include "sixDOFODE.H"
#include "tolerancesSwitch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class geometricSixDOF Declaration
\*---------------------------------------------------------------------------*/
class geometricSixDOF
:
public sixDOFODE
{
// Private data
// Initial body state variables
//- Displacement relative to spring equilibrium
dimensionedVector Xrel_;
//- Velocity of mass centroid
dimensionedVector U_;
//- Average velocity of mass centroid (evaluated at midstep)
dimensionedVector Uaverage_;
//- Rotation tensor
tensor rotation_;
//- Rotational increment tensor
tensor rotIncrement_;
//- Rotational velocity about mass centroid
dimensionedVector omega_;
//- Average rotational velocity in relative coordinate system
// (evaluated at midstep)
dimensionedVector omegaAverage_;
//- ODE coefficients
scalarField coeffs_;
// Private Member Functions
//- Disallow default bitwise copy construct
geometricSixDOF(const geometricSixDOF&);
//- Disallow default bitwise assignment
void operator=(const geometricSixDOF&);
// Variables in relative coordinate system (solved for)
//- Return acceleration in relative coordinate system
// given current values of relative displacement, velocity,
// orientation (rotation tensor) and time
dimensionedVector A
(
const dimensionedVector& xR,
const dimensionedVector& uR,
const tensor& R,
const scalar t
) const;
//- Return rotational acceleration in relative coordinate system
// given current values for relative rotational velocity,
// orientation (rotation tensor) and time
dimensionedVector OmegaDot
(
const tensor& R,
const dimensionedVector& omega,
const scalar t
) const;
//- Return the Euler part of moment equation
dimensionedVector E
(
const dimensionedVector& omega
) const;
//- Exponential map used to calculate increment of the rotation
// tensor
tensor expMap(const vector& rotInc) const;
//- Differential of the expontential map used to calculate the time
// derivative of rotation increment vector
vector dexpMap(const vector& rotInc, const vector& omega) const;
protected:
// Protected Member Functions
// Non-access control for setting state variables
//- Set ODE parameters from another ODE
virtual void setState(const sixDOFODE&);
public:
// Run-time type information
TypeName("geometricSixDOF");
// Static data members
//- Rotational increment tensor tolerance. Used in expMap member
// function in case the rotation is negligibly small
static const debug::tolerancesSwitch rotIncTensorTol_;
//- Rotational increment rate of change tolerance. Used in dexpMap
// member function in case the rotation rate is negligibly small
static const debug::tolerancesSwitch rotIncRateTol_;
// Constructors
//- Construct from dictionary
geometricSixDOF(const IOobject& io);
//- Construct geometricSixDOF object, changing name
geometricSixDOF
(
const word& name,
const geometricSixDOF& gsd
);
//- Return a clone, changing the name
virtual autoPtr<sixDOFODE> clone(const word& name) const;
// Destructor
virtual ~geometricSixDOF();
// Member Functions
// Virtual interface for 6DOF motion state
// Variables in relative coordinate system
//- Return displacement in translated coordinate system
// relative to spring equilibrium
virtual const dimensionedVector& Xrel() const;
//- Return rotational velocity in relative coordinate system
virtual const dimensionedVector& omega() const;
// Displacement and velocity in the absolute coordinate system
//- Return position of origin in absolute coordinate system
virtual dimensionedVector X() const;
//- Return velocity of origin
virtual const dimensionedVector& U() const;
//- Return average velocity of origin (evaluated at midstep)
virtual const dimensionedVector& Uaverage() const;
// Average motion per time-step
//- Return average rotational velocity in relative coordinate
// system (evaluated at midstep)
virtual const dimensionedVector& omegaAverage() const;
// Rotations
//- Return rotation tensor to relative coordinate system
virtual tensor toRelative() const;
//- Return rotation tensor to absolute coordinate system
virtual tensor toAbsolute() const;
//- Return transformation tensor between new and previous
// rotation
virtual const tensor& rotIncrementTensor() const;
// ODE parameters
//- Return number of equations
virtual label nEqns() const
{
return 12;
}
//- Return access to coefficients
virtual scalarField& coeffs()
{
return coeffs_;
}
//- Return reference to coefficients
virtual const scalarField& coeffs() const
{
return coeffs_;
}
//- Evaluate derivatives
virtual void derivatives
(
const scalar x,
const scalarField& y,
scalarField& dydx
) const;
//- Evaluate Jacobian
virtual void jacobian
(
const scalar x,
const scalarField& y,
scalarField& dfdx,
scalarSquareMatrix& dfdy
) const
{
notImplemented
(
"geometricSixDOF::jacobian\n"
"(\n"
" const scalar x,\n"
" const scalarField& y,\n"
" scalarField& dfdx,\n"
" scalarSquareMatrix& dfdy,\n"
") const"
);
}
//- Update ODE after the solution, advancing by delta
virtual void update(const scalar delta);
// Write controls
//- WriteData member function required by regIOobject
virtual bool writeData(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

Some files were not shown because too many files have changed in this diff Show more