diff --git a/ThirdParty/AllMake.stage2 b/ThirdParty/AllMake.stage2 index 44613e46c..b8d04fdeb 100755 --- a/ThirdParty/AllMake.stage2 +++ b/ThirdParty/AllMake.stage2 @@ -40,8 +40,8 @@ # Requirements: # 1: Your foam-extend environment must be properly initialized # 2: AllMake.stage1 if you are overriding your system compiler -# 3: The file etc/prefs.sh should be used for setting the variables enabling -# the compilation of the various packages +# 3: The file etc/prefs.sh should be used for setting the variables +# enabling the compilation of the various packages # # Author: # Martin Beaudoin, Hydro-Quebec, (2015) @@ -104,6 +104,13 @@ then ( rpm_make -p openmpi-1.8.8 -s openmpi-1.8.8.spec -u http://www.open-mpi.org/software/ompi/v1.8/downloads/openmpi-1.8.8.tar.gz \ -f --define '_configureAdditionalArgs "$WM_THIRD_PARTY_USE_OPENMPI_188_ConfigureAdditionalArgs"') } + # mvipich2-2.2 + # + [ ! -z "$WM_THIRD_PARTY_USE_MVAPICH2_22" ] && { + echo "Building mvapich2 2.2" + ( rpm_make -p mvapich2-2.2 -s mvapich2-2.2.spec -u http://mvapich.cse.ohio-state.edu/download/mvapich/mv2/mvapich2-2.2.tar.gz \ + -f --define '_configureAdditionalArgs "$WM_THIRD_PARTY_USE_MVAPICH2_22_ConfigureAdditionalArgs"') + } else echo "Using system installed OpenMPI" echo "" diff --git a/ThirdParty/rpmBuild/SPECS/mvapich2-2.2.spec b/ThirdParty/rpmBuild/SPECS/mvapich2-2.2.spec new file mode 100644 index 000000000..fed9345c4 --- /dev/null +++ b/ThirdParty/rpmBuild/SPECS/mvapich2-2.2.spec @@ -0,0 +1,268 @@ +#------------------------------------------------------------------------------ +# ========= | +# \\ / F ield | foam-extend: Open Source CFD +# \\ / O peration | +# \\ / A nd | For copyright notice see file Copyright +# \\/ M anipulation | +#------------------------------------------------------------------------------ +# 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 . +# +# Script +# RPM spec file for mvapich2-2.2 +# +# Description +# RPM spec file for creating a relocatable RPM +# +# Author: +# Martin Beaudoin, Hydro-Quebec, (2015) +# +#------------------------------------------------------------------------------ + +# We grab the value of WM_THIRD_PARTY and WM_OPTIONS from the environment variable +%{expand:%%define _WM_THIRD_PARTY_DIR %(echo $WM_THIRD_PARTY_DIR)} +%{expand:%%define _WM_OPTIONS %(echo $WM_OPTIONS)} + +# Disable the generation of debuginfo packages +%define debug_package %{nil} + +# The topdir needs to point to the $WM_THIRD_PARTY/rpmbuild directory +%define _topdir %{_WM_THIRD_PARTY_DIR}/rpmBuild +%define _tmppath %{_topdir}/tmp + +# Will install the package directly $WM_THIRD_PARTY_DIR +# Some comments about package relocation: +# By using this prefix for the Prefix: parameter in this file, you will make this +# package relocatable. +# +# This is fine, as long as your software is itself relocatable. +# +# Simply take note that libraries built with libtool are not relocatable because the +# prefix we specify will be hard-coded in the library .la files. +# Ref: http://sourceware.org/autobook/autobook/autobook_80.html +# +# In that case, if you ever change the value of the $WM_THIRD_PARTY_DIR, you will +# not be able to reutilize this RPM, even though it is relocatable. You will need to +# regenerate the RPM. +# +%define _prefix %{_WM_THIRD_PARTY_DIR} + +%define name mvapich2 +%define release %{_WM_OPTIONS} +%define version 2.2 + +%define buildroot %{_topdir}/BUILD/%{name}-%{version}-root + +BuildRoot: %{buildroot} +Summary: mvapich2 +License: BSD +Name: %{name} +Version: %{version} +Release: %{release} +URL: http://mvapich.cse.ohio-state.edu/download/mvapich/mv2/ +Source: %url/%{name}-%{version}.tar.gz +Prefix: %{_prefix} +Group: Development/Tools + + +%define _installPrefix %{_prefix}/packages/%{name}-%{version}/platforms/%{_WM_OPTIONS} + +#-------------------------------------------------------------------------- +# +# Here, we define default compiling options for mvapich2 +# +# One can override the option on the commande line : --define='MACRO EXPR' +# +%{!?_configureAdditionalArgs: %define _configureAdditionalArgs Undefined} + + +%description +%{summary} + +%prep +%setup -q + + +%build + # export WM settings in a form that GNU configure recognizes + [ -n "$WM_CC" ] && export CC="$WM_CC" + [ -n "$WM_CXX" ] && export CXX="$WM_CXX" + [ -n "$WM_CFLAGS" ] && export CFLAGS="$WM_CFLAGS" + [ -n "$WM_CXXFLAGS" ] && export CXXFLAGS="$WM_CXXFLAGS" + [ -n "$WM_LDFLAGS" ] && export LDFLAGS="$WM_LDFLAGS" + + set +x + echo "" + echo "Compilation options:" + echo " _configureAdditionalArgs : %{_configureAdditionalArgs}" + echo "" + set -x + + unset mpiWith + # Enable GridEngine if it appears to be in use + # If you don't want any integration with SGE, simply unset the SGE + # environment variable + if [ -n "$SGE_ROOT" ] + then + mpiWith="$mpiWith --with-sge" + else + mpiWith="$mpiWith --without-sge" + mpiWith="$mpiWith --enable-mca-no-build=ras-gridengine,pls-gridengine" + fi + + # Infiniband support + # Adjust according to your local paths + # if [ -d /usr/local/ofed -a -d /usr/local/ofed/lib64 ] + # then + # mpiWith="$mpiWith --with-openib=/usr/local/ofed" + # mpiWith="$mpiWith --with-openib-libdir=/usr/local/ofed/lib64" + # fi + + # NOTE. Option --disable-mcast disables Infiniband when libs are not present + + ./configure \ + --prefix=%{_installPrefix} \ + --exec_prefix=%{_installPrefix} \ + --enable-mpirun-prefix-by-default \ + --enable-orterun-prefix-by-default \ + --enable-shared \ + --enable-mpi-cxx \ + --disable-static \ + --disable-mpi-f77 \ + --disable-mpi-f90 \ + --without-slurm \ + --enable-mpi-profile $mpiWith \ + --disable-vt \ + echo %{?_configureAdditionalArgs}` + + [ -z "$WM_NCOMPPROCS" ] && WM_NCOMPPROCS=1 + make -j $WM_NCOMPPROCS + +%install + make install DESTDIR=$RPM_BUILD_ROOT + + # Creation of foam-extend specific .csh and .sh files" + + echo "" + echo "Generating foam-extend specific .csh and .sh files for the package %{name}-%{version}" + echo "" + # + # Generate package specific .sh file for foam-extend + # +mkdir -p $RPM_BUILD_ROOT/%{_installPrefix}/etc +cat << DOT_SH_EOF > $RPM_BUILD_ROOT/%{_installPrefix}/etc/%{name}-%{version}.sh +# Load %{name}-%{version} libraries and binaries +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +export MVAPICH2_DIR=\$WM_THIRD_PARTY_DIR/packages/%{name}-%{version}/platforms/\$WM_OPTIONS +export MVAPICH2_BIN_DIR=\$MVAPICH2_DIR/bin +export MVAPICH2_LIB_DIR=\$MVAPICH2_DIR/lib + +# Enable access to the package runtime applications and libraries +[ -d \$MVAPICH2_BIN_DIR ] && _foamAddPath \$MVAPICH2_BIN_DIR +[ -d \$MVAPICH2_LIB_DIR ] && _foamAddLib \$MVAPICH2_LIB_DIR + +export MPI_HOME=\$MVAPICH2_DIR +export MPI_ARCH_PATH=\$MPI_HOME +export OPAL_PREFIX=\$MPI_ARCH_PATH + +# We initialize the rest of the environment using mpicc --showme: +# This does not work with MVAPICH2 +#export MVAPICH2_INCLUDE_DIR="\`mpicc --showme:incdirs\`" +#export MVAPICH2_COMPILE_FLAGS="\`mpicc --showme:compile\`" +#export MVAPICH2_LINK_FLAGS="\`mpicc --showme:link\`" + +# Set the foam-extend compilation flags +export PINC=\$MVAPICH2_COMPILE_FLAGS +export PLIBS=\$MVAPICH2_LINK_FLAGS + + +if [ "\$FOAM_VERBOSE" -a "\$PS1" ] +then + echo " Environment variables defined for OpenMPI:" + echo " MVAPICH2_BIN_DIR : \$MVAPICH2_BIN_DIR" + echo " MVAPICH2_LIB_DIR : \$MVAPICH2_LIB_DIR" + echo " MVAPICH2_INCLUDE_DIR : \$MVAPICH2_INCLUDE_DIR" + echo " MVAPICH2_COMPILE_FLAGS : \$MVAPICH2_COMPILE_FLAGS" + echo " MVAPICH2_LINK_FLAGS : \$MVAPICH2_LINK_FLAGS" + echo "" + echo " MPI_HOME : \$MPI_HOME" + echo " MPI_ARCH_PATH : \$MPI_ARCH_PATH" + echo " OPAL_PREFIX : \$OPAL_PREFIX" + echo " PINC : \$PINC" + echo " PLIBS : \$PLIBS" +fi +DOT_SH_EOF + + # + # Generate package specific .csh file for foam-extend + # +cat << DOT_CSH_EOF > $RPM_BUILD_ROOT/%{_installPrefix}/etc/%{name}-%{version}.csh +# Load %{name}-%{version} libraries and binaries +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +setenv MVAPICH2_DIR \$WM_THIRD_PARTY_DIR/packages/%{name}-%{version}/platforms/\$WM_OPTIONS +setenv MVAPICH2_BIN_DIR \$MVAPICH2_DIR/bin +setenv MVAPICH2_LIB_DIR \$MVAPICH2_DIR/lib + +# Enable access to the package runtime applications and libraries +if ( -e \$MVAPICH2_BIN_DIR ) then + _foamAddPath \$MVAPICH2_BIN_DIR +endif + +if ( -e \$MVAPICH2_LIB_DIR ) then + _foamAddLib \$MVAPICH2_LIB_DIR +endif + +setenv MPI_HOME \$MVAPICH2_DIR +setenv MPI_ARCH_PATH \$MPI_HOME +setenv OPAL_PREFIX \$MPI_ARCH_PATH + +# We initialize the rest of the environment using mpicc --showme: +setenv MVAPICH2_INCLUDE_DIR "\`mpicc --showme:incdirs\`" +setenv MVAPICH2_COMPILE_FLAGS "\`mpicc --showme:compile\`" +setenv MVAPICH2_LINK_FLAGS "\`mpicc --showme:link\`" + +# Set the foam-extend compilation flags +setenv PINC "\$MVAPICH2_COMPILE_FLAGS" +setenv PLIBS "\$MVAPICH2_LINK_FLAGS" + + +if (\$?FOAM_VERBOSE && \$?prompt) then + echo " Environment variables defined for OpenMPI:" + echo " MVAPICH2_BIN_DIR : \$MVAPICH2_BIN_DIR" + echo " MVAPICH2_LIB_DIR : \$MVAPICH2_LIB_DIR" + echo " MVAPICH2_INCLUDE_DIR : \$MVAPICH2_INCLUDE_DIR" + echo " MVAPICH2_COMPILE_FLAGS : \$MVAPICH2_COMPILE_FLAGS" + echo " MVAPICH2_LINK_FLAGS : \$MVAPICH2_LINK_FLAGS" + echo "" + echo " MPI_HOME : \$MPI_HOME" + echo " MPI_ARCH_PATH : \$MPI_ARCH_PATH" + echo " OPAL_PREFIX : \$OPAL_PREFIX" + echo " PINC : \$PINC" + echo " PLIBS : \$PLIBS" +endif +DOT_CSH_EOF + + #finally, generate a .tgz file for systems where using rpm for installing packages + # as a non-root user might be a problem. + (mkdir -p %{_topdir}/TGZS/%{_target_cpu}; cd $RPM_BUILD_ROOT/%{_prefix}; tar -zcvf %{_topdir}/TGZS/%{_target_cpu}/%{name}-%{version}.tgz packages/%{name}-%{version}) + + +%clean +rm -rf %{buildroot} + +%files +%defattr(-,root,root) +%{_installPrefix} diff --git a/ThirdParty/rpmBuild/SPECS/openmpi-1.8.8.spec b/ThirdParty/rpmBuild/SPECS/openmpi-1.8.8.spec index 78530a615..913d2e1cc 100644 --- a/ThirdParty/rpmBuild/SPECS/openmpi-1.8.8.spec +++ b/ThirdParty/rpmBuild/SPECS/openmpi-1.8.8.spec @@ -138,7 +138,7 @@ Group: Development/Tools --enable-orterun-prefix-by-default \ --enable-shared \ --enable-mpi-cxx \ - --disable-static \ + --disable-static \ --disable-mpi-f77 \ --disable-mpi-f90 \ --without-slurm \ diff --git a/ThirdParty/tools/makeThirdPartyFunctionsForRPM b/ThirdParty/tools/makeThirdPartyFunctionsForRPM index b256f3944..f2308f4e0 100755 --- a/ThirdParty/tools/makeThirdPartyFunctionsForRPM +++ b/ThirdParty/tools/makeThirdPartyFunctionsForRPM @@ -82,6 +82,7 @@ rpm_make() echo "Package name : $_PACKAGE" echo "Package URL : $_PACKAGE_URL" echo "RPM spec file name: $_SPECFILE" + echo "RPM file name : $_RPMFILENAME" echo "Additional flags : $_ADDITIONALFLAGS" # Shift options diff --git a/applications/solvers/coupled/MRFPorousFoam/UEqn.H b/applications/solvers/coupled/MRFPorousFoam/UEqn.H index 80e70f3a5..c40bb8ef2 100644 --- a/applications/solvers/coupled/MRFPorousFoam/UEqn.H +++ b/applications/solvers/coupled/MRFPorousFoam/UEqn.H @@ -5,6 +5,9 @@ fvc::div(phi) ); + // Update boundary velocity for consistency with the flux + mrfZones.correctBoundaryVelocity(U); + // Momentum equation fvVectorMatrix UEqn ( @@ -67,28 +70,38 @@ const scalarField& V = mesh.V().field(); - // Note: this insertion should only happen in porous cell zones - // A rewrite of the porousZones class interface is needed. + // Note: insertion should only happen in porous cell zones // HJ, 14/Mar/2016 - forAll (TUIn, cellI) + + register label cellI; + + forAll (pZones, pZoneI) { - const scalar& cellV = V[cellI]; + const labelList& curZoneCells = pZones[pZoneI].zone(); - const tensor& cellTU = TUIn[cellI]; + // Loop over all cells in the zone + forAll (curZoneCells, zcI) + { + cellI = curZoneCells[zcI]; - CoeffField::squareType& cellDD = DD[cellI]; + const scalar& cellV = V[cellI]; - cellDD(0, 0) += cellV*cellTU.xx(); - cellDD(0, 1) += cellV*cellTU.xy(); - cellDD(0, 2) += cellV*cellTU.xz(); + const tensor& cellTU = TUIn[cellI]; - cellDD(1, 0) += cellV*cellTU.yx(); - cellDD(1, 1) += cellV*cellTU.yy(); - cellDD(2, 2) += cellV*cellTU.yz(); + CoeffField::squareType& cellDD = DD[cellI]; - cellDD(2, 0) += cellV*cellTU.zx(); - cellDD(2, 1) += cellV*cellTU.zy(); - cellDD(2, 2) += cellV*cellTU.zz(); + cellDD(0, 0) += cellV*cellTU.xx(); + cellDD(0, 1) += cellV*cellTU.xy(); + cellDD(0, 2) += cellV*cellTU.xz(); + + cellDD(1, 0) += cellV*cellTU.yx(); + cellDD(1, 1) += cellV*cellTU.yy(); + cellDD(2, 2) += cellV*cellTU.yz(); + + cellDD(2, 0) += cellV*cellTU.zx(); + cellDD(2, 1) += cellV*cellTU.zy(); + cellDD(2, 2) += cellV*cellTU.zz(); + } } } } diff --git a/applications/solvers/incompressible/MRFSimpleFoam/MRFSimpleFoam.C b/applications/solvers/incompressible/MRFSimpleFoam/MRFSimpleFoam.C index a17092cad..996fb9631 100644 --- a/applications/solvers/incompressible/MRFSimpleFoam/MRFSimpleFoam.C +++ b/applications/solvers/incompressible/MRFSimpleFoam/MRFSimpleFoam.C @@ -61,54 +61,14 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector { - // Momentum predictor - - tmp UEqn - ( - fvm::div(phi, U) - + turbulence->divDevReff() - ); - mrfZones.addCoriolis(UEqn()); - - UEqn().relax(); - - solve(UEqn() == -fvc::grad(p)); - - volScalarField rAU = 1.0/UEqn().A(); - U = rAU*UEqn().H(); - UEqn.clear(); - - phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf(); - mrfZones.relativeFlux(phi); - adjustPhi(phi, U, p); - - // Non-orthogonal pressure corrector loop - while (simple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - fvm::laplacian(rAU, p) == fvc::div(phi) - ); - - pEqn.setReference(pRefCell, pRefValue); - pEqn.solve(); - - if (simple.finalNonOrthogonalIter()) - { - phi -= pEqn.flux(); - } - } - -# include "continuityErrs.H" - - // Explicitly relax pressure for momentum corrector - p.relax(); - - // Momentum corrector - U -= rAU*fvc::grad(p); - U.correctBoundaryConditions(); +# include "UEqn.H" +# include "pEqn.H" } + // Calculate relative velocity + Urel == U; + mrfZones.relativeVelocity(Urel); + turbulence->correct(); runTime.write(); diff --git a/applications/solvers/incompressible/MRFSimpleFoam/UEqn.H b/applications/solvers/incompressible/MRFSimpleFoam/UEqn.H new file mode 100644 index 000000000..db0cdfb73 --- /dev/null +++ b/applications/solvers/incompressible/MRFSimpleFoam/UEqn.H @@ -0,0 +1,19 @@ + // Solve the momentum equation + + tmp HUEqn + ( + fvm::div(phi, U) + + turbulence->divDevReff() + ); + mrfZones.addCoriolis(HUEqn()); + + // Get under-relaxation factor + const scalar UUrf = mesh.solutionDict().equationRelaxationFactor(U.name()); + + // Momentum solution + solve + ( + relax(HUEqn(), UUrf) + == + -fvc::grad(p) + ); diff --git a/applications/solvers/incompressible/MRFSimpleFoam/createFields.H b/applications/solvers/incompressible/MRFSimpleFoam/createFields.H index 8ae490c37..e43f40963 100644 --- a/applications/solvers/incompressible/MRFSimpleFoam/createFields.H +++ b/applications/solvers/incompressible/MRFSimpleFoam/createFields.H @@ -44,3 +44,19 @@ MRFZones mrfZones(mesh); mrfZones.correctBoundaryVelocity(U); + + // Create Urel as a permanent field to make it available for on-the-fly + // post-processing operations + volVectorField Urel + ( + IOobject + ( + "Urel", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + U + ); + mrfZones.relativeVelocity(Urel); diff --git a/applications/solvers/incompressible/MRFSimpleFoam/pEqn.H b/applications/solvers/incompressible/MRFSimpleFoam/pEqn.H new file mode 100644 index 000000000..1752fc663 --- /dev/null +++ b/applications/solvers/incompressible/MRFSimpleFoam/pEqn.H @@ -0,0 +1,50 @@ +{ + p.boundaryField().updateCoeffs(); + + // Prepare clean 1/Ap without contribution from under-relaxation + // HJ, 26/Oct/2015 + volScalarField rAU + ( + "(1|A(U))", + 1/HUEqn().A() + ); + + // Store velocity under-relaxation point before using U for + // the flux precursor + U.storePrevIter(); + + U = rAU*HUEqn().H(); + HUEqn.clear(); + + phi = fvc::interpolate(U) & mesh.Sf(); + mrfZones.relativeFlux(phi); + adjustPhi(phi, U, p); + + // Non-orthogonal pressure corrector loop + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(rAU, p) == fvc::div(phi) + ); + + pEqn.setReference(pRefCell, pRefValue); + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi -= pEqn.flux(); + } + } + +# include "continuityErrs.H" + + // Explicitly relax pressure for momentum corrector + p.relax(); + + // Momentum corrector + // Note: since under-relaxation does not change aU, H/a in U can be + // re-used. HJ, 22/Jan/2016 + U = UUrf*(U - rAU*fvc::grad(p)) + (1 - UUrf)*U.prevIter(); + U.correctBoundaryConditions(); +} diff --git a/applications/solvers/incompressible/simpleFoam/pEqn.H b/applications/solvers/incompressible/simpleFoam/pEqn.H index 9a1c0a4ef..aede8a3ba 100644 --- a/applications/solvers/incompressible/simpleFoam/pEqn.H +++ b/applications/solvers/incompressible/simpleFoam/pEqn.H @@ -1,3 +1,4 @@ +{ p.boundaryField().updateCoeffs(); // Prepare clean 1/Ap without contribution from under-relaxation @@ -14,6 +15,7 @@ U = rUA*HUEqn().H(); HUEqn.clear(); + phi = fvc::interpolate(U) & mesh.Sf(); adjustPhi(phi, U, p); @@ -22,7 +24,7 @@ { fvScalarMatrix pEqn ( - fvm::laplacian(rUA, p) == fvc::div(phi) + fvm::laplacian(rAU, p) == fvc::div(phi) ); pEqn.setReference(pRefCell, pRefValue); @@ -43,6 +45,6 @@ // Momentum corrector // Note: since under-relaxation does not change aU, H/a in U can be // re-used. HJ, 22/Jan/2016 - U = UUrf*(U - rUA*fvc::grad(p)) + (1 - UUrf)*U.prevIter(); + U = UUrf*(U - rAU*fvc::grad(p)) + (1 - UUrf)*U.prevIter(); U.correctBoundaryConditions(); - +} diff --git a/bin/plotForces.py b/bin/plotForces.py index ff59ca41a..a7d779848 100755 --- a/bin/plotForces.py +++ b/bin/plotForces.py @@ -53,6 +53,19 @@ for i in range(1,len(t)): my[i] += mvy[i] mz[i] += mvz[i] +# write clean data file +outForces=open('forces.dat','w') + +for data in zip(t,fx,fy,fz): + outForces.write(' '.join([str(d) for d in data])+'\n') +outForces.close() + +outMoments=open('moments.dat','w') + +for data in zip(t,mx,my,mz): + outMoments.write(' '.join([str(d) for d in data])+'\n') +outMoments.close() + # plot forces import pylab pylab.xlabel('iteration') diff --git a/bin/plotResidual.py b/bin/plotResidual.py index 22149b34e..5122812d4 100755 --- a/bin/plotResidual.py +++ b/bin/plotResidual.py @@ -62,6 +62,7 @@ for line in lines: tepsilon.append(iepsilon) epsilon.append(float(matchepsilon.group(2))) +# write clean data file outfile=open('residual.dat','w') if iomega > 0: diff --git a/etc/settings.csh b/etc/settings.csh index fd84b9d46..7ef3dce6e 100644 --- a/etc/settings.csh +++ b/etc/settings.csh @@ -300,12 +300,17 @@ case SYSTEMOPENMPI: breaksw case MVAPICH2: - set mpi_version=mvapich2-2.2 - if ($?WM_THIRD_PARTY_USE_MVAPICH2_22 != 0 && -d $WM_THIRD_PARTY_DIR/packages/mvapich2-2.2/platforms/$WM_OPTIONS ) then - if ($?FOAM_VERBOSE && $?prompt) then - echo "Using mvapich2-2.2 from the ThirdParty package: $WM_THIRD_PARTY_DIR/packages/$mpi_version" + set mpi_version=mvapich2 + + if ($?MVAPICH2_BIN_DIR != 0) then + if (-d "${MVAPICH2_BIN_DIR}" ) then + _foamAddPath $MVAPICH2_BIN_DIR endif - _foamSource $WM_THIRD_PARTY_DIR/packages/$mpi_version/platforms/$WM_OPTIONS/etc/$mpi_version.csh + else + set mpicc_cmd=`which mpicc` + setenv MVAPICH2_BIN_DIR `dirname $mpicc_cmd` + unset mpicc_cmd + endif setenv MPI_HOME `dirname $MVAPICH2_BIN_DIR` setenv MPI_ARCH_PATH $MPI_HOME @@ -571,8 +576,8 @@ endif # zoltan # ~~~~~ -if ( $?ZOLTAN_SYSTEM == 0 && $?WM_THIRD_PARTY_USE_ZOLTAN_36 != 0 && -e "$WM_THIRD_PARTY_DIR"/packages/zoltan-3.6/platforms/$WM_OPTIONS ) then - _foamSource $WM_THIRD_PARTY_DIR/packages/zoltan-3.6/platforms/$WM_OPTIONS/etc/zoltan-3.6.csh +if ( $?ZOLTAN_SYSTEM == 0 && $?WM_THIRD_PARTY_USE_ZOLTAN_35 != 0 && -e "$WM_THIRD_PARTY_DIR"/packages/zoltan-3.5/platforms/$WM_OPTIONS ) then + _foamSource $WM_THIRD_PARTY_DIR/packages/zoltan-3.5/platforms/$WM_OPTIONS/etc/zoltan-3.5.csh endif # Python diff --git a/etc/settings.sh b/etc/settings.sh index eb59fa877..32c636fc8 100755 --- a/etc/settings.sh +++ b/etc/settings.sh @@ -375,22 +375,16 @@ SYSTEMOPENMPI) unset mpi_version ;; + MVAPICH2) - mpi_version=mvapich2 - - if [ -n "${MVAPICH2_BIN_DIR}" ] && [ -d "${MVAPICH2_BIN_DIR}" ] - then - _foamAddPath $MVAPICH2_BIN_DIR - else - MVAPICH2_BIN_DIR=$(dirname `which mpicc`) - fi - - if which mpicc >/dev/null - then - mpicc -v 2>/dev/null | grep -q "mpicc for MVAPICH2" || - echo "Warning: `which mpicc` does not appear to be for MVAPICH2" - else - echo "Warning: mpicc not available" + mpi_version=mvapich2-2.2 + if [ ! -z $WM_THIRD_PARTY_USE_MVAPICH2_20 ] && [ -e $WM_THIRD_PARTY_DIR/packages/mvapich2-2.2/platforms/$WM_OPTIONS ] + then + if [ "$FOAM_VERBOSE" -a "$PS1" ] + then + echo "Using mvapich2 from the ThirdParty package: $WM_THIRD_PARTY_DIR/packages/$mpi_version" + fi + _foamSource $WM_THIRD_PARTY_DIR/packages/$mpi_version/platforms/$WM_OPTIONS/etc/$mpi_version.sh fi export MPI_HOME=`dirname $MVAPICH2_BIN_DIR` diff --git a/src/finiteArea/Make/files b/src/finiteArea/Make/files index 6cd6d6288..9e8ffce86 100644 --- a/src/finiteArea/Make/files +++ b/src/finiteArea/Make/files @@ -48,6 +48,8 @@ derivedFaPatchFields = $(faPatchFields)/derived $(derivedFaPatchFields)/fixedValueOutflow/fixedValueOutflowFaPatchFields.C $(derivedFaPatchFields)/inletOutlet/inletOutletFaPatchFields.C $(derivedFaPatchFields)/slip/slipFaPatchFields.C +$(derivedFaPatchFields)/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.C +$(derivedFaPatchFields)/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.C faePatchFields = fields/faePatchFields $(faePatchFields)/faePatchField/faePatchFields.C @@ -92,6 +94,7 @@ $(ddtSchemes)/backwardFaDdtScheme/backwardFaDdtSchemes.C $(ddtSchemes)/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.C divSchemes = finiteArea/divSchemes +finiteArea/fam/vectorFamDiv.C $(divSchemes)/faDivScheme/faDivSchemes.C $(divSchemes)/gaussFaDivScheme/gaussFaDivSchemes.C diff --git a/src/finiteArea/faMesh/faMesh.C b/src/finiteArea/faMesh/faMesh.C index 481ccaf35..1c1ab187b 100644 --- a/src/finiteArea/faMesh/faMesh.C +++ b/src/finiteArea/faMesh/faMesh.C @@ -1047,6 +1047,18 @@ void Foam::faMesh::addFaPatches(const List& p) } +Foam::label Foam::faMesh::comm() const +{ + return comm_; +} + + +Foam::label& Foam::faMesh::comm() +{ + return comm_; +} + + const Foam::objectRegistry& Foam::faMesh::thisDb() const { return mesh().thisDb(); diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H index 8f340701b..89f81abec 100644 --- a/src/finiteArea/faMesh/faMesh.H +++ b/src/finiteArea/faMesh/faMesh.H @@ -110,6 +110,12 @@ class faMesh mutable label nFaces_; + // Communication support + + //- Communicator used for parallel communication + label comm_; + + // Demand-driven data //- Primitive patch @@ -287,9 +293,8 @@ public: ); - // Destructor - - virtual ~faMesh(); + //- Destructor + virtual ~faMesh(); // Member Functions @@ -369,6 +374,15 @@ public: } + // Communication support + + //- Return communicator used for parallel communication + label comm() const; + + //- Return communicator used for parallel communication + label& comm(); + + // Access //- Return reference to the mesh database diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C index 050f763d1..08c4809bb 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C +++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C @@ -55,6 +55,18 @@ processorFaPatch::~processorFaPatch() // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // +Foam::label Foam::processorFaPatch::comm() const +{ + return boundaryMesh().mesh().comm(); +} + + +int Foam::processorFaPatch::tag() const +{ + return Pstream::msgType(); +} + + void processorFaPatch::makeNonGlobalPatchPoints() const { // If it is not runing parallel or there are no global points diff --git a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H index 6939342aa..0f0c96d81 100644 --- a/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H +++ b/src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H @@ -37,7 +37,6 @@ SourceFiles #include "coupledFaPatch.H" #include "processorLduInterface.H" -// #include "processorPolyPatch.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -55,7 +54,10 @@ class processorFaPatch { // Private data + //- My processro number int myProcNo_; + + //- Neighbour processor number int neighbProcNo_; //- Processor-neighbbour patch edge centres @@ -75,6 +77,7 @@ class processorFaPatch // non-global, i.e. present in this processor patch mutable labelList* nonGlobalPatchPointsPtr_; + protected: // Protected Member functions @@ -88,27 +91,27 @@ protected: //- Find non-globa patch points void makeNonGlobalPatchPoints() const; -protected: - // Protected Member functions + // Geometry functions - //- Initialise the calculation of the patch geometry - void initGeometry(); + //- Initialise the calculation of the patch geometry + void initGeometry(); - //- Calculate the patch geometry - void calcGeometry(); + //- Calculate the patch geometry + void calcGeometry(); - //- Initialise the patches for moving points - void initMovePoints(const pointField&); + //- Initialise the patches for moving points + void initMovePoints(const pointField&); - //- Correct patches after moving points - void movePoints(const pointField&); + //- Correct patches after moving points + void movePoints(const pointField&); - //- Initialise the update of the patch topology - virtual void initUpdateMesh(); + //- Initialise the update of the patch topology + virtual void initUpdateMesh(); + + //- Update of the patch topology + virtual void updateMesh(); - //- Update of the patch topology - virtual void updateMesh(); public: @@ -160,9 +163,9 @@ public: nonGlobalPatchPointsPtr_(NULL) {} - // Destructor - virtual ~processorFaPatch(); + //- Destructor + virtual ~processorFaPatch(); // Member functions @@ -192,6 +195,16 @@ public: } } + + // Communications support + + //- Return communicator used for communication + virtual label comm() const; + + //- Return message tag to use for communication + virtual int tag() const; + + //- Return face transformation tensor virtual const tensorField& forwardT() const { diff --git a/src/finiteArea/fields/faPatchFields/basic/transform/transformFaPatchFields.C b/src/finiteArea/fields/faPatchFields/basic/transform/transformFaPatchFields.C index 2c808f7bd..24e9fbfc8 100644 --- a/src/finiteArea/fields/faPatchFields/basic/transform/transformFaPatchFields.C +++ b/src/finiteArea/fields/faPatchFields/basic/transform/transformFaPatchFields.C @@ -28,7 +28,6 @@ Description #include "faPatchFields.H" #include "transformFaPatchFields.H" #include "addToRunTimeSelectionTable.H" -// #include "volFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -37,10 +36,6 @@ namespace Foam // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -// defineNamedTemplateTypeNameAndDebug(transformFaPatchScalarField, 0); -// defineNamedTemplateTypeNameAndDebug(transformFaPatchVectorField, 0); -// defineNamedTemplateTypeNameAndDebug(transformFaPatchTensorField, 0); - makeFaPatchFieldsTypeName(transform); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C b/src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C index 4e155a0a6..3537b4c27 100644 --- a/src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C +++ b/src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchField.C @@ -176,7 +176,7 @@ void processorFaPatchField::initEvaluate { if (Pstream::parRun()) { - procPatch_.compressedSend(commsType, this->patchInternalField()()); + procPatch_.send(commsType, this->patchInternalField()()); } } @@ -189,7 +189,7 @@ void processorFaPatchField::evaluate { if (Pstream::parRun()) { - procPatch_.compressedReceive(commsType, *this); + procPatch_.receive(commsType, *this); if (doTransform()) { @@ -218,7 +218,7 @@ void processorFaPatchField::initInterfaceMatrixUpdate const bool switchToLhs ) const { - procPatch_.compressedSend + procPatch_.send ( commsType, this->patch().patchInternalField(psiInternal)() @@ -240,7 +240,7 @@ void processorFaPatchField::updateInterfaceMatrix { scalarField pnf ( - procPatch_.compressedReceive(commsType, this->size())() + procPatch_.receive(commsType, this->size())() ); // Transform according to the transformation tensor diff --git a/src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C b/src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C index 11be0915b..fe5b4dc27 100644 --- a/src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C +++ b/src/finiteArea/fields/faPatchFields/constraint/processor/processorFaPatchScalarField.C @@ -52,7 +52,7 @@ void processorFaPatchField::initInterfaceMatrixUpdate const bool switchToLhs ) const { - procPatch_.compressedSend + procPatch_.send ( commsType, patch().patchInternalField(psiInternal)() @@ -74,7 +74,7 @@ void processorFaPatchField::updateInterfaceMatrix { scalarField pnf ( - procPatch_.compressedReceive(commsType, this->size())() + procPatch_.receive(commsType, this->size())() ); const unallocLabelList& edgeFaces = patch().edgeFaces(); diff --git a/src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchFields.C b/src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchFields.C index b67ab97c2..be4183c9a 100644 --- a/src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchFields.C +++ b/src/finiteArea/fields/faPatchFields/constraint/wedge/wedgeFaPatchFields.C @@ -28,7 +28,6 @@ Description #include "faPatchFields.H" #include "wedgeFaPatchFields.H" #include "addToRunTimeSelectionTable.H" -// #include "volFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -37,10 +36,6 @@ namespace Foam // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -// defineNamedTemplateTypeNameAndDebug(wedgeFaPatchScalarField, 0); -// defineNamedTemplateTypeNameAndDebug(wedgeFaPatchVectorField, 0); -// defineNamedTemplateTypeNameAndDebug(wedgeFaPatchTensorField, 0); - makeFaPatchFields(wedge); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteArea/fields/faPatchFields/derived/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.C b/src/finiteArea/fields/faPatchFields/derived/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.C new file mode 100644 index 000000000..bc0a93108 --- /dev/null +++ b/src/finiteArea/fields/faPatchFields/derived/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.C @@ -0,0 +1,155 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "edgeNormalFixedValueFaPatchVectorField.H" +#include "addToRunTimeSelectionTable.H" +#include "areaFields.H" +#include "faPatchFieldMapper.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::edgeNormalFixedValueFaPatchVectorField:: +edgeNormalFixedValueFaPatchVectorField +( + const faPatch& p, + const DimensionedField& iF +) +: + fixedValueFaPatchVectorField(p, iF), + refValue_(p.size(), 0) +{} + + +Foam::edgeNormalFixedValueFaPatchVectorField:: +edgeNormalFixedValueFaPatchVectorField +( + const edgeNormalFixedValueFaPatchVectorField& ptf, + const faPatch& p, + const DimensionedField& iF, + const faPatchFieldMapper& mapper +) +: + fixedValueFaPatchVectorField(ptf, p, iF, mapper), + refValue_(ptf.refValue_, mapper) +{} + + +Foam::edgeNormalFixedValueFaPatchVectorField:: +edgeNormalFixedValueFaPatchVectorField +( + const faPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFaPatchVectorField(p, iF, dict), + refValue_("refValue", dict, p.size()) +{} + + +Foam::edgeNormalFixedValueFaPatchVectorField:: +edgeNormalFixedValueFaPatchVectorField +( + const edgeNormalFixedValueFaPatchVectorField& pivpvf +) +: + fixedValueFaPatchVectorField(pivpvf), + refValue_(pivpvf.refValue_) +{} + + +Foam::edgeNormalFixedValueFaPatchVectorField:: +edgeNormalFixedValueFaPatchVectorField +( + const edgeNormalFixedValueFaPatchVectorField& pivpvf, + const DimensionedField& iF +) +: + fixedValueFaPatchVectorField(pivpvf, iF), + refValue_(pivpvf.refValue_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::edgeNormalFixedValueFaPatchVectorField::autoMap +( + const faPatchFieldMapper& m +) +{ + fixedValueFaPatchVectorField::autoMap(m); + refValue_.autoMap(m); +} + + +void Foam::edgeNormalFixedValueFaPatchVectorField::rmap +( + const faPatchVectorField& ptf, + const labelList& addr +) +{ + fixedValueFaPatchVectorField::rmap(ptf, addr); + + const edgeNormalFixedValueFaPatchVectorField& tiptf = + refCast(ptf); + + refValue_.rmap(tiptf.refValue_, addr); +} + + +void Foam::edgeNormalFixedValueFaPatchVectorField::updateCoeffs() +{ + if (updated()) + { + return; + } + + // Bug fix: update for moving mesh. HJ, 15/Oct/2010 + operator==(refValue_*patch().edgeNormals()); +} + + +void Foam::edgeNormalFixedValueFaPatchVectorField::write(Ostream& os) const +{ + fixedValueFaPatchVectorField::write(os); + refValue_.writeEntry("refValue", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +makeFaPatchTypeField +( + faPatchVectorField, + edgeNormalFixedValueFaPatchVectorField +); + +} // End namespace Foam + + +// ************************************************************************* // diff --git a/src/finiteArea/fields/faPatchFields/derived/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.H b/src/finiteArea/fields/faPatchFields/derived/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.H new file mode 100644 index 000000000..0cf215d73 --- /dev/null +++ b/src/finiteArea/fields/faPatchFields/derived/edgeNormalFixedValue/edgeNormalFixedValueFaPatchVectorField.H @@ -0,0 +1,173 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::edgeNormalFixedValueFaPatchVectorField + +Description + Edge normal fixed value vector field finite area boundary condition + + Describes a surface normal vector boundary condition by its magnitude. + Note: The value is positive for outward-pointing vectors + +SourceFiles + edgeNormalFixedValueFaPatchVectorField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef edgeNormalFixedValueFaPatchVectorField_H +#define edgeNormalFixedValueFaPatchVectorField_H + +#include "faPatchFields.H" +#include "fixedValueFaPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class edgeNormalFixedValueFaPatch Declaration +\*---------------------------------------------------------------------------*/ + +class edgeNormalFixedValueFaPatchVectorField +: + public fixedValueFaPatchVectorField +{ + // Private data + + //- Surface-normal velocity value + scalarField refValue_; + + +public: + + //- Runtime type information + TypeName("edgeNormalFixedValue"); + + + // Constructors + + //- Construct from patch and internal field + edgeNormalFixedValueFaPatchVectorField + ( + const faPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + edgeNormalFixedValueFaPatchVectorField + ( + const faPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given + // edgeNormalFixedValueFaPatchVectorField + // onto a new patch + edgeNormalFixedValueFaPatchVectorField + ( + const edgeNormalFixedValueFaPatchVectorField&, + const faPatch&, + const DimensionedField&, + const faPatchFieldMapper& + ); + + //- Construct as copy + edgeNormalFixedValueFaPatchVectorField + ( + const edgeNormalFixedValueFaPatchVectorField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new edgeNormalFixedValueFaPatchVectorField(*this) + ); + } + + //- Construct as copy setting internal field reference + edgeNormalFixedValueFaPatchVectorField + ( + const edgeNormalFixedValueFaPatchVectorField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new edgeNormalFixedValueFaPatchVectorField + ( + *this, + iF + ) + ); + } + + + + // Member functions + + // Mapping functions + + //- Map (and resize as needed) from self given a mapping object + virtual void autoMap + ( + const faPatchFieldMapper& + ); + + //- Reverse map the given faPatchField onto this faPatchField + virtual void rmap + ( + const faPatchVectorField&, + const labelList& + ); + + + // Member functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchField.C b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchField.C new file mode 100644 index 000000000..cb2b83c26 --- /dev/null +++ b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchField.C @@ -0,0 +1,137 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "timeVaryingUniformFixedValueFaPatchField.H" +#include "foamTime.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::timeVaryingUniformFixedValueFaPatchField:: +timeVaryingUniformFixedValueFaPatchField +( + const faPatch& p, + const DimensionedField& iF +) +: + fixedValueFaPatchField(p, iF), + timeSeries_() +{} + + +template +Foam::timeVaryingUniformFixedValueFaPatchField:: +timeVaryingUniformFixedValueFaPatchField +( + const faPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFaPatchField(p, iF), + timeSeries_(dict) +{ + if (dict.found("value")) + { + faPatchField::operator==(Field("value", dict, p.size())); + } + else + { + updateCoeffs(); + } +} + + +template +Foam::timeVaryingUniformFixedValueFaPatchField:: +timeVaryingUniformFixedValueFaPatchField +( + const timeVaryingUniformFixedValueFaPatchField& ptf, + const faPatch& p, + const DimensionedField& iF, + const faPatchFieldMapper& mapper +) +: + fixedValueFaPatchField(ptf, p, iF, mapper), + timeSeries_(ptf.timeSeries_) +{} + + +template +Foam::timeVaryingUniformFixedValueFaPatchField:: +timeVaryingUniformFixedValueFaPatchField +( + const timeVaryingUniformFixedValueFaPatchField& ptf +) +: + fixedValueFaPatchField(ptf), + timeSeries_(ptf.timeSeries_) +{} + + +template +Foam::timeVaryingUniformFixedValueFaPatchField:: +timeVaryingUniformFixedValueFaPatchField +( + const timeVaryingUniformFixedValueFaPatchField& ptf, + const DimensionedField& iF +) +: + fixedValueFaPatchField(ptf, iF), + timeSeries_(ptf.timeSeries_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::timeVaryingUniformFixedValueFaPatchField::updateCoeffs() +{ + if (this->updated()) + { + return; + } + + faPatchField::operator== + ( + timeSeries_(this->db().time().timeOutputValue()) + ); + fixedValueFaPatchField::updateCoeffs(); +} + + +template +void Foam::timeVaryingUniformFixedValueFaPatchField::write +( + Ostream& os +) const +{ + faPatchField::write(os); + timeSeries_.write(os); + this->writeEntry("value", os); +} + + +// ************************************************************************* // diff --git a/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchField.H b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchField.H new file mode 100644 index 000000000..b59f386c4 --- /dev/null +++ b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchField.H @@ -0,0 +1,182 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::timeVaryingUniformFixedValueFaPatchField + +Description + A time-varying form of a uniform fixed value finite area + boundary condition. + + Example of the boundary condition specification: + @verbatim + inlet + { + type timeVaryingUniformFixedValue; + fileName "$FOAM_CASE/time-series"; + outOfBounds clamp; // (error|warn|clamp|repeat) + } + @endverbatim + +Note + This class is derived directly from a fixedValue patch rather than from + a uniformFixedValue patch. + +See Also + Foam::interpolationTable and Foam::fixedValueFaPatchField + +SourceFiles + timeVaryingUniformFixedValueFaPatchField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef timeVaryingUniformFixedValueFaPatchField_H +#define timeVaryingUniformFixedValueFaPatchField_H + +#include "fixedValueFaPatchField.H" +#include "interpolationTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class timeVaryingUniformFixedValueFaPatch Declaration +\*---------------------------------------------------------------------------*/ + +template +class timeVaryingUniformFixedValueFaPatchField +: + public fixedValueFaPatchField +{ + // Private data + + //- The time series being used, including the bounding treatment + interpolationTable timeSeries_; + + +public: + + //- Runtime type information + TypeName("timeVaryingUniformFixedValue"); + + + // Constructors + + //- Construct from patch and internal field + timeVaryingUniformFixedValueFaPatchField + ( + const faPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + timeVaryingUniformFixedValueFaPatchField + ( + const faPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given patch field onto a new patch + timeVaryingUniformFixedValueFaPatchField + ( + const timeVaryingUniformFixedValueFaPatchField&, + const faPatch&, + const DimensionedField&, + const faPatchFieldMapper& + ); + + //- Construct as copy + timeVaryingUniformFixedValueFaPatchField + ( + const timeVaryingUniformFixedValueFaPatchField& + ); + + //- Construct and return a clone + virtual tmp > clone() const + { + return tmp > + ( + new timeVaryingUniformFixedValueFaPatchField(*this) + ); + } + + //- Construct as copy setting internal field reference + timeVaryingUniformFixedValueFaPatchField + ( + const timeVaryingUniformFixedValueFaPatchField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp > clone + ( + const DimensionedField& iF + ) const + { + return tmp > + ( + new timeVaryingUniformFixedValueFaPatchField(*this, iF) + ); + } + + + // Member functions + + // Access + + //- Return the time series used + const interpolationTable& timeSeries() const + { + return timeSeries_; + } + + + // Evaluation functions + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "timeVaryingUniformFixedValueFaPatchField.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/foam/db/IOstreams/Pstreams/PstreamsPrint.C b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.C similarity index 76% rename from src/foam/db/IOstreams/Pstreams/PstreamsPrint.C rename to src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.C index ac0fdf58c..0610f7e76 100644 --- a/src/foam/db/IOstreams/Pstreams/PstreamsPrint.C +++ b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.C @@ -21,28 +21,23 @@ License You should have received a copy of the GNU General Public License along with foam-extend. If not, see . -Description - Prints out a description of the streams - \*---------------------------------------------------------------------------*/ -#include "IPstream.H" -#include "OPstream.H" +#include "timeVaryingUniformFixedValueFaPatchFields.H" +#include "addToRunTimeSelectionTable.H" +#include "areaFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -void Foam::IPstream::print(Ostream& os) const +namespace Foam { - os << "Reading from processor " << fromProcNo_ - << " to processor " << myProcNo() << Foam::endl; -} +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -void Foam::OPstream::print(Ostream& os) const -{ - os << "Writing from processor " << toProcNo_ - << " to processor " << myProcNo() << Foam::endl; -} +makeFaPatchFields(timeVaryingUniformFixedValue); +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam // ************************************************************************* // diff --git a/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.H b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.H new file mode 100644 index 000000000..71f70dd63 --- /dev/null +++ b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFields.H @@ -0,0 +1,49 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#ifndef timeVaryingUniformFixedValueFaPatchFields_H +#define timeVaryingUniformFixedValueFaPatchFields_H + +#include "timeVaryingUniformFixedValueFaPatchField.H" +#include "fieldTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +makeFaPatchTypeFieldTypedefs(timeVaryingUniformFixedValue) + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFieldsFwd.H b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFieldsFwd.H new file mode 100644 index 000000000..2a8360f27 --- /dev/null +++ b/src/finiteArea/fields/faPatchFields/derived/timeVaryingUniformFixedValue/timeVaryingUniformFixedValueFaPatchFieldsFwd.H @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#ifndef timeVaryingUniformFixedValueFaPatchFieldsFwd_H +#define timeVaryingUniformFixedValueFaPatchFieldsFwd_H + +#include "faPatchField.H" +#include "fieldTypes.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template class timeVaryingUniformFixedValueFaPatchField; + +makeFaPatchTypeFieldTypedefs(timeVaryingUniformFixedValue); + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.H b/src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.H index dc8fefadd..018b276ac 100644 --- a/src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.H +++ b/src/finiteArea/fields/faPatchFields/faPatchField/faPatchField.H @@ -468,60 +468,87 @@ public: #endif -#define makeFaPatchTypeFieldTypeName(typePatchTypeField) \ - \ -defineNamedTemplateTypeNameAndDebug(typePatchTypeField, 0); - -#define makeFaPatchFieldsTypeName(typePatchField) \ - \ -makeFaPatchTypeFieldTypeName(typePatchField##FaPatchScalarField); \ -makeFaPatchTypeFieldTypeName(typePatchField##FaPatchVectorField); \ -makeFaPatchTypeFieldTypeName(typePatchField##FaPatchSphericalTensorField); \ -makeFaPatchTypeFieldTypeName(typePatchField##FaPatchSymmTensorField); \ -makeFaPatchTypeFieldTypeName(typePatchField##FaPatchTensorField); - -#define makeFaPatchTypeField(PatchTypeField, typePatchTypeField) \ - \ -defineNamedTemplateTypeNameAndDebug(typePatchTypeField, 0); \ - \ -addToRunTimeSelectionTable \ -( \ - PatchTypeField, typePatchTypeField, patch \ -); \ - \ -addToRunTimeSelectionTable \ -( \ - PatchTypeField, \ - typePatchTypeField, \ - patchMapper \ -); \ - \ -addToRunTimeSelectionTable \ -( \ - PatchTypeField, typePatchTypeField, dictionary \ +#define addToFaPatchFieldRunTimeSelection(PatchTypeField, typePatchTypeField) \ + \ +addToRunTimeSelectionTable \ +( \ + PatchTypeField, typePatchTypeField, patch \ +); \ + \ +addToRunTimeSelectionTable \ +( \ + PatchTypeField, \ + typePatchTypeField, \ + patchMapper \ +); \ + \ +addToRunTimeSelectionTable \ +( \ + PatchTypeField, typePatchTypeField, dictionary \ ); -#define makeFaPatchFields(type) \ - \ -makeFaPatchTypeField(faPatchScalarField, type##FaPatchScalarField); \ -makeFaPatchTypeField(faPatchVectorField, type##FaPatchVectorField); \ -makeFaPatchTypeField \ -( \ - faPatchSphericalTensorField, \ - type##FaPatchSphericalTensorField \ -); \ -makeFaPatchTypeField(faPatchSymmTensorField, type##FaPatchSymmTensorField);\ -makeFaPatchTypeField(faPatchTensorField, type##FaPatchTensorField); +#define makeFaPatchTypeFieldTypeName(typePatchTypeField) \ + \ +defineNamedTemplateTypeNameAndDebug(typePatchTypeField, 0); -#define makeFaPatchTypeFieldTypedefs(type) \ - \ -typedef type##FaPatchField type##FaPatchScalarField; \ -typedef type##FaPatchField type##FaPatchVectorField; \ -typedef type##FaPatchField \ - type##FaPatchSphericalTensorField; \ -typedef type##FaPatchField type##FaPatchSymmTensorField; \ +#define makeFaPatchFieldsTypeName(typePatchField) \ + \ +makeFaPatchTypeFieldTypeName(typePatchField##FaPatchScalarField); \ +makeFaPatchTypeFieldTypeName(typePatchField##FaPatchVectorField); \ +makeFaPatchTypeFieldTypeName(typePatchField##FaPatchSphericalTensorField); \ +makeFaPatchTypeFieldTypeName(typePatchField##FaPatchSymmTensorField); \ +makeFaPatchTypeFieldTypeName(typePatchField##FaPatchTensorField); + + +#define makeFaPatchTypeField(PatchTypeField, typePatchTypeField) \ + \ +defineTypeNameAndDebug(typePatchTypeField, 0); \ + \ +addToFaPatchFieldRunTimeSelection \ +( \ + PatchTypeField, typePatchTypeField \ +); + +#define makeTemplateFaPatchTypeField(PatchTypeField, typePatchTypeField) \ + \ +defineNamedTemplateTypeNameAndDebug(typePatchTypeField, 0); \ + \ +addToFaPatchFieldRunTimeSelection \ +( \ + PatchTypeField, typePatchTypeField \ +); + + +#define makeFaPatchFields(type) \ + \ +makeTemplateFaPatchTypeField(faPatchScalarField, type##FaPatchScalarField); \ +makeTemplateFaPatchTypeField(faPatchVectorField, type##FaPatchVectorField); \ +makeTemplateFaPatchTypeField \ +( \ + faPatchSphericalTensorField, \ + type##FaPatchSphericalTensorField \ +); \ +makeTemplateFaPatchTypeField \ +( \ + faPatchSymmTensorField, \ + type##FaPatchSymmTensorField \ +); \ +makeTemplateFaPatchTypeField \ +( \ + faPatchTensorField, \ + type##FaPatchTensorField \ +); + + +#define makeFaPatchTypeFieldTypedefs(type) \ + \ +typedef type##FaPatchField type##FaPatchScalarField; \ +typedef type##FaPatchField type##FaPatchVectorField; \ +typedef type##FaPatchField \ + type##FaPatchSphericalTensorField; \ +typedef type##FaPatchField type##FaPatchSymmTensorField; \ typedef type##FaPatchField type##FaPatchTensorField; diff --git a/src/finiteArea/finiteArea/divSchemes/gaussFaDivScheme/gaussFaDivScheme.C b/src/finiteArea/finiteArea/divSchemes/gaussFaDivScheme/gaussFaDivScheme.C index c4d43aefc..802b543e0 100644 --- a/src/finiteArea/finiteArea/divSchemes/gaussFaDivScheme/gaussFaDivScheme.C +++ b/src/finiteArea/finiteArea/divSchemes/gaussFaDivScheme/gaussFaDivScheme.C @@ -62,7 +62,7 @@ gaussDivScheme::facDiv ( this->mesh_.Le() & this->tinterpScheme_().interpolate(vf) ) - - this->mesh_.faceCurvatures()*(this->mesh_.faceAreaNormals()&vf) + // Removed for consistencty. Matthias Rauter, 6/Dec/2016 ); tDiv().rename("div(" + vf.name() + ')'); diff --git a/src/finiteArea/finiteArea/fac/fac.H b/src/finiteArea/finiteArea/fac/fac.H index ce5e382b9..2b17490bc 100644 --- a/src/finiteArea/finiteArea/fac/fac.H +++ b/src/finiteArea/finiteArea/fac/fac.H @@ -43,6 +43,8 @@ Description #include "facAverage.H" #include "facLnGrad.H" #include "facDdt.H" +#include "facNGrad.H" +#include "facNDiv.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteArea/finiteArea/fac/facDiv.C b/src/finiteArea/finiteArea/fac/facDiv.C index cc72a13d9..e5f04f9d9 100644 --- a/src/finiteArea/finiteArea/fac/facDiv.C +++ b/src/finiteArea/finiteArea/fac/facDiv.C @@ -28,6 +28,7 @@ License #include "facEdgeIntegrate.H" #include "faDivScheme.H" #include "faConvectionScheme.H" +#include "transformField.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -48,7 +49,17 @@ div const GeometricField& ssf ) { - return fac::edgeIntegrate(ssf); + const areaVectorField& n = ssf.mesh().faceAreaNormals(); + + tmp > tDiv = + fac::edgeIntegrate(ssf); + + GeometricField& Div = tDiv(); + + Div.internalField() = transform(tensor::I - sqr(n), Div.internalField()); + Div.correctBoundaryConditions(); + + return tDiv; } @@ -79,10 +90,34 @@ div const word& name ) { - return fa::divScheme::New + const areaVectorField& n = vf.mesh().faceAreaNormals(); + + tmp + < + GeometricField + < + typename innerProduct::type, + faPatchField, + areaMesh + > + > tDiv ( - vf.mesh(), vf.mesh().schemesDict().divScheme(name) - )().facDiv(vf); + fa::divScheme::New + ( + vf.mesh(), vf.mesh().schemesDict().divScheme(name) + )().facDiv(vf) + ); + GeometricField + < + typename innerProduct::type, + faPatchField, + areaMesh + >& Div = tDiv(); + + Div.internalField() = transform(tensor::I - sqr(n), Div.internalField()); + Div.correctBoundaryConditions(); + + return tDiv; } @@ -159,12 +194,24 @@ div const word& name ) { - return fa::convectionScheme::New + const areaVectorField& n = vf.mesh().faceAreaNormals(); + + tmp > tDiv ( - vf.mesh(), - flux, - vf.mesh().schemesDict().divScheme(name) - )().facDiv(flux, vf); + fa::convectionScheme::New + ( + vf.mesh(), + flux, + vf.mesh().schemesDict().divScheme(name) + )().facDiv(flux, vf) + ); + GeometricField& Div = tDiv(); + + Div.internalField() = transform(tensor::I - sqr(n), Div.internalField()); + Div.correctBoundaryConditions(); + + return tDiv; + } diff --git a/src/finiteArea/finiteArea/fac/facGrad.C b/src/finiteArea/finiteArea/fac/facGrad.C index e7c3b5e88..568f4c2b6 100644 --- a/src/finiteArea/finiteArea/fac/facGrad.C +++ b/src/finiteArea/finiteArea/fac/facGrad.C @@ -55,7 +55,18 @@ grad const GeometricField& ssf ) { - return fac::edgeIntegrate(ssf.mesh().Sf() * ssf); + const areaVectorField &n = ssf.mesh().faceAreaNormals(); + typedef typename outerProduct::type GradType; + + tmp > tgGrad = + fac::edgeIntegrate(ssf.mesh().Sf()*ssf); + + GeometricField& gGrad = tgGrad(); + + gGrad -= (gGrad & n)*n; + gGrad.correctBoundaryConditions(); + + return tgGrad; } template @@ -95,11 +106,22 @@ grad const word& name ) { - return fa::gradScheme::New - ( - vf.mesh(), - vf.mesh().schemesDict().gradScheme(name) - )().grad(vf); + const areaVectorField &n = vf.mesh().faceAreaNormals(); + typedef typename outerProduct::type GradType; + + tmp > tgGrad = + fa::gradScheme::New + ( + vf.mesh(), + vf.mesh().schemesDict().gradScheme(name) + )().grad(vf); + + GeometricField& gGrad = tgGrad(); + + gGrad -= (gGrad & n)*n; + gGrad.correctBoundaryConditions(); + + return tgGrad; } diff --git a/src/finiteArea/finiteArea/fac/facNDiv.C b/src/finiteArea/finiteArea/fac/facNDiv.C new file mode 100644 index 000000000..6a1953961 --- /dev/null +++ b/src/finiteArea/finiteArea/fac/facNDiv.C @@ -0,0 +1,333 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "facNDiv.H" +#include "faMesh.H" +#include "facEdgeIntegrate.H" +#include "faDivScheme.H" +#include "faConvectionScheme.H" +#include "transformField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fac +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +tmp > +ndiv +( + const GeometricField& ssf +) +{ + const areaVectorField &n = ssf.mesh().faceAreaNormals(); + + tmp > v = fac::edgeIntegrate(ssf); + + //v.internalField() = transform(n*n, v.internalField()); + v.internalField() = (v.internalField()&n)*n; + + return v; +} + + +template +tmp > +ndiv +( + const tmp >& tssf +) +{ + tmp > Div(fac::ndiv(tssf())); + tssf.clear(); + return Div; +} + + +template +tmp +< + GeometricField + < + typename innerProduct::type, faPatchField, areaMesh + > +> +ndiv +( + const GeometricField& vf, + const word& name +) +{ + const areaVectorField &n = vf.mesh().faceAreaNormals(); + + tmp > tDiv + ( + fa::divScheme::New + ( + vf.mesh(), vf.mesh().schemesDict().divScheme(name) + )().facDiv(vf) + ); + + GeometricField& Div = tDiv(); + + //Div.internalField() = transform(n*n, Div.internalField()); + Div.internalField() = (Div.internalField()&n)*n; + return tDiv; +} + + +template +tmp +< + GeometricField + < + typename innerProduct::type, faPatchField, areaMesh + > +> +ndiv +( + const tmp >& tvvf, + const word& name +) +{ + typedef typename innerProduct::type DivType; + tmp > Div + ( + fac::ndiv(tvvf(), name) + ); + tvvf.clear(); + return Div; +} + +template +tmp +< + GeometricField + < + typename innerProduct::type, faPatchField, areaMesh + > +> +ndiv +( + const GeometricField& vf +) +{ + return fac::ndiv(vf, "div("+vf.name()+')'); +} + + +template +tmp +< + GeometricField + < + typename innerProduct::type, faPatchField, areaMesh + > +> +ndiv +( + const tmp >& tvvf +) +{ + typedef typename innerProduct::type DivType; + tmp > Div + ( + fac::ndiv(tvvf()) + ); + + tvvf.clear(); + return Div; +} + + +template +tmp > +ndiv +( + const edgeScalarField& flux, + const GeometricField& vf, + const word& name +) +{ + const areaVectorField &n = vf.mesh().faceAreaNormals(); + + tmp > tDiv + ( + fa::convectionScheme::New + ( + vf.mesh(), + flux, + vf.mesh().schemesDict().divScheme(name) + )().facDiv(flux, vf) + ); + + GeometricField& Div = tDiv(); + + //Div.internalField() = transform(n*n, Div.internalField()); + Div.internalField() = (Div.internalField()&n)*n; + + return tDiv; + +} + + +template +tmp > +ndiv +( + const tmp& tflux, + const GeometricField& vf, + const word& name +) +{ + tmp > Div + ( + fac::ndiv(tflux(), vf, name) + ); + tflux.clear(); + return Div; +} + + +template +tmp > +ndiv +( + const edgeScalarField& flux, + const tmp >& tvf, + const word& name +) +{ + tmp > Div + ( + fac::ndiv(flux, tvf(), name) + ); + tvf.clear(); + return Div; +} + + +template +tmp > +ndiv +( + const tmp& tflux, + const tmp >& tvf, + const word& name +) +{ + tmp > Div + ( + fac::ndiv(tflux(), tvf(), name) + ); + tflux.clear(); + tvf.clear(); + return Div; +} + + +template +tmp > +ndiv +( + const edgeScalarField& flux, + const GeometricField& vf +) +{ + return fac::ndiv + ( + flux, vf, "div("+flux.name()+','+vf.name()+')' + ); +} + + +template +tmp > +ndiv +( + const tmp& tflux, + const GeometricField& vf +) +{ + tmp > Div + ( + fac::ndiv(tflux(), vf) + ); + tflux.clear(); + return Div; +} + + +template +tmp > +ndiv +( + const edgeScalarField& flux, + const tmp >& tvf +) +{ + tmp > Div + ( + fac::ndiv(flux, tvf()) + ); + tvf.clear(); + return Div; +} + + +template +tmp > +ndiv +( + const tmp& tflux, + const tmp >& tvf +) +{ + tmp > Div + ( + fac::ndiv(tflux(), tvf()) + ); + tflux.clear(); + tvf.clear(); + return Div; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fac + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteArea/finiteArea/fac/facNDiv.H b/src/finiteArea/finiteArea/fac/facNDiv.H new file mode 100644 index 000000000..7aab1da79 --- /dev/null +++ b/src/finiteArea/finiteArea/fac/facNDiv.H @@ -0,0 +1,187 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Namespace + fac + +Description + Calculate the divergence of the given field. + +SourceFiles + facDiv.C + +\*---------------------------------------------------------------------------*/ + + +#ifndef facNDiv_H +#define facNDiv_H + +#include "areaFieldsFwd.H" +#include "edgeFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Namespace fac functions Declaration +\*---------------------------------------------------------------------------*/ + +namespace fac +{ + template + tmp > ndiv + ( + const GeometricField& + ); + + template + tmp > ndiv + ( + const tmp >& + ); + + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + > ndiv + ( + const GeometricField&, + const word& name + ); + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + > ndiv + ( + const tmp >&, + const word& name + ); + + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + > ndiv + ( + const GeometricField& + ); + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + > ndiv + ( + const tmp >& + ); + + + template + tmp > ndiv + ( + const edgeScalarField&, + const GeometricField&, + const word& name + ); + + template + tmp > ndiv + ( + const tmp&, + const GeometricField&, + const word& name + ); + + template + tmp > ndiv + ( + const edgeScalarField&, + const tmp >&, + const word& name + ); + + template + tmp > ndiv + ( + const tmp&, + const tmp >&, + const word& name + ); + + + template + tmp > ndiv + ( + const edgeScalarField&, + const GeometricField& + ); + + template + tmp > ndiv + ( + const tmp&, + const GeometricField& + ); + + template + tmp > ndiv + ( + const edgeScalarField&, + const tmp >& + ); + + template + tmp > ndiv + ( + const tmp&, + const tmp >& + ); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "facNDiv.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/finiteArea/fac/facNGrad.C b/src/finiteArea/finiteArea/fac/facNGrad.C new file mode 100644 index 000000000..9d4d96070 --- /dev/null +++ b/src/finiteArea/finiteArea/fac/facNGrad.C @@ -0,0 +1,202 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "facNGrad.H" +#include "areaFields.H" +#include "edgeFields.H" +#include "facEdgeIntegrate.H" +#include "faMesh.H" +#include "faGradScheme.H" +#include "transformField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fac +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +tmp +< + GeometricField + < + typename outerProduct::type, faPatchField, areaMesh + > +> +ngrad +( + const GeometricField& ssf +) +{ + const areaVectorField &n = ssf.mesh().faceAreaNormals(); + typedef typename outerProduct::type GradType; + + tmp > tgGrad = fac::edgeIntegrate(ssf.mesh().Sf() * ssf); + + GeometricField& grad = tgGrad(); + + grad = (grad&n)*n; + + return tgGrad; +} + +template +tmp +< + GeometricField + < + typename outerProduct::type, faPatchField, areaMesh + > +> +ngrad +( + const tmp >& tssf +) +{ + typedef typename outerProduct::type GradType; + tmp > Grad + ( + fac::ngrad(tssf()) + ); + tssf.clear(); + return Grad; +} + + +template +tmp +< + GeometricField + < + typename outerProduct::type, faPatchField, areaMesh + > +> +ngrad +( + const GeometricField& vf, + const word& name +) +{ + const areaVectorField &n = vf.mesh().faceAreaNormals(); + typedef typename outerProduct::type GradType; + + tmp > tgGrad = fa::gradScheme::New + ( + vf.mesh(), + vf.mesh().schemesDict().gradScheme(name) + )().grad(vf); + + GeometricField& grad = tgGrad(); + + grad = (grad&n)*n; + + return tgGrad; +} + + +template +tmp +< + GeometricField + < + typename outerProduct::type, faPatchField, areaMesh + > +> +ngrad +( + const tmp >& tvf, + const word& name +) +{ + tmp + < + GeometricField + < + typename outerProduct::type, faPatchField, areaMesh + > + > tGrad + ( + fac::ngrad(tvf(), name) + ); + tvf.clear(); + return tGrad; +} + + +template +tmp +< + GeometricField + < + typename outerProduct::type, faPatchField, areaMesh + > +> +ngrad +( + const GeometricField& vf +) +{ + return fac::ngrad(vf, "grad(" + vf.name() + ')'); +} + + +template +tmp +< + GeometricField + < + typename outerProduct::type, faPatchField, areaMesh + > +> +ngrad +( + const tmp >& tvf +) +{ + typedef typename outerProduct::type GradType; + tmp > Grad + ( + fac::ngrad(tvf()) + ); + tvf.clear(); + return Grad; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fac + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteArea/finiteArea/fac/facNGrad.H b/src/finiteArea/finiteArea/fac/facNGrad.H new file mode 100644 index 000000000..0375202bd --- /dev/null +++ b/src/finiteArea/finiteArea/fac/facNGrad.H @@ -0,0 +1,131 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Namespace + fac + +Description + Calculate the gradient normal to the surface of the given field. + +SourceFiles + facGrad.C + +\*---------------------------------------------------------------------------*/ + + +#ifndef facNGrad_H +#define facNGrad_H + +#include "areaFieldsFwd.H" +#include "edgeFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Namespace fac functions Declaration +\*---------------------------------------------------------------------------*/ + +namespace fac +{ + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + > ngrad + ( + const GeometricField& + ); + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + > ngrad + ( + const tmp >& + ); + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + >ngrad + ( + const GeometricField&, + const word& name + ); + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + >ngrad + ( + const tmp >&, + const word& name + ); + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + >ngrad + ( + const GeometricField& + ); + + template + tmp + < + GeometricField + ::type, faPatchField, areaMesh> + >ngrad + ( + const tmp >& + ); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "facNGrad.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/finiteArea/fam/fam.H b/src/finiteArea/finiteArea/fam/fam.H index 65b1c6905..5e9f7ef33 100644 --- a/src/finiteArea/finiteArea/fam/fam.H +++ b/src/finiteArea/finiteArea/fam/fam.H @@ -43,6 +43,7 @@ SourceFiles #include "famDiv.H" #include "famLaplacian.H" #include "famSup.H" +#include "famNDiv.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/finiteArea/finiteArea/fam/famDiv.C b/src/finiteArea/finiteArea/fam/famDiv.C index b382b60a5..2296b224b 100644 --- a/src/finiteArea/finiteArea/fam/famDiv.C +++ b/src/finiteArea/finiteArea/fam/famDiv.C @@ -49,14 +49,36 @@ div const word& name ) { - return fa::convectionScheme::New + const areaVectorField& n = vf.mesh().faceAreaNormals(); + + tmp > tM ( - vf.mesh(), - flux, - vf.mesh().schemesDict().divScheme(name) - )().famDiv(flux, vf); + fa::convectionScheme::New + ( + vf.mesh(), + flux, + vf.mesh().schemesDict().divScheme(name) + )().famDiv(flux, vf) + ); + faMatrix& M = tM(); + + GeometricField v + ( + fa::convectionScheme::New + ( + vf.mesh(), + flux, + vf.mesh().schemesDict().divScheme(name) + )().facDiv(flux, vf) + ); + + //HJ Check if the product is from left or right. HJ, 6/Dec/2016 + M -= (v & n)*n; + + return tM; } + template tmp > div diff --git a/src/finiteArea/finiteArea/fam/famDiv.H b/src/finiteArea/finiteArea/fam/famDiv.H index 33b92f49d..e8c4f2bb2 100644 --- a/src/finiteArea/finiteArea/fam/famDiv.H +++ b/src/finiteArea/finiteArea/fam/famDiv.H @@ -29,6 +29,7 @@ Description SourceFiles famDiv.C + vectorFamDiv.C \*---------------------------------------------------------------------------*/ @@ -89,6 +90,9 @@ namespace fam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Template specialisation +#include "vectorFamDiv.H" + #ifdef NoRepository # include "famDiv.C" #endif diff --git a/src/finiteArea/finiteArea/fam/famNDiv.C b/src/finiteArea/finiteArea/fam/famNDiv.C new file mode 100644 index 000000000..33c55f3fa --- /dev/null +++ b/src/finiteArea/finiteArea/fam/famNDiv.C @@ -0,0 +1,108 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "famNDiv.H" +#include "faMesh.H" +#include "faMatrix.H" +#include "faConvectionScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +tmp > +ndiv +( + const edgeScalarField& flux, + GeometricField& vf, + const word& name +) +{ + return fa::convectionScheme::New + ( + vf.mesh(), + flux, + vf.mesh().schemesDict().divScheme(name) + )().famDiv(flux, vf);//TODO calculate normal +} + +template +tmp > +ndiv +( + const tmp& tflux, + GeometricField& vf, + const word& name +) +{ + tmp > Div(fam::ndiv(tflux(), vf, name)); + tflux.clear(); + return Div; +} + + +template +tmp > +ndiv +( + const edgeScalarField& flux, + GeometricField& vf +) +{ + return fam::ndiv(flux, vf, "div("+flux.name()+','+vf.name()+')'); +} + +template +tmp > +ndiv +( + const tmp& tflux, + GeometricField& vf +) +{ + tmp > Div(fam::ndiv(tflux(), vf)); + tflux.clear(); + return Div; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteArea/finiteArea/fam/famNDiv.H b/src/finiteArea/finiteArea/fam/famNDiv.H new file mode 100644 index 000000000..0abf49e54 --- /dev/null +++ b/src/finiteArea/finiteArea/fam/famNDiv.H @@ -0,0 +1,100 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 3.2 + \\ / 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 . + +Namespace + fam + +Description + Calculate the matrix for the divergence of the given field and flux. + +SourceFiles + famDiv.C + +\*---------------------------------------------------------------------------*/ + +#ifndef famNDiv_H +#define famNDiv_H + +#include "areaFieldsFwd.H" +#include "edgeFieldsFwd.H" +#include "edgeInterpolationScheme.H" +#include "faMatrix.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Namespace fam functions Declaration +\*---------------------------------------------------------------------------*/ + +namespace fam +{ + template + tmp > ndiv + ( + const edgeScalarField&, + GeometricField&, + const word& name + ); + + template + tmp > ndiv + ( + const tmp&, + GeometricField&, + const word& name + ); + + template + tmp > ndiv + ( + const edgeScalarField&, + GeometricField& + ); + + template + tmp > ndiv + ( + const tmp&, + GeometricField& + ); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "famNDiv.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/finiteArea/fam/vectorFamDiv.C b/src/finiteArea/finiteArea/fam/vectorFamDiv.C new file mode 100644 index 000000000..c3b941a3e --- /dev/null +++ b/src/finiteArea/finiteArea/fam/vectorFamDiv.C @@ -0,0 +1,69 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "vectorFamDiv.H" +#include "faMesh.H" +#include "faMatrix.H" +#include "faConvectionScheme.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace fam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<> +tmp > +div +( + const edgeScalarField& flux, + const GeometricField& vf, + const word& name +) +{ + return fa::convectionScheme::New + ( + vf.mesh(), + flux, + vf.mesh().schemesDict().divScheme(name) + )().famDiv(flux, vf); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace fam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/finiteArea/finiteArea/fam/vectorFamDiv.H b/src/finiteArea/finiteArea/fam/vectorFamDiv.H new file mode 100644 index 000000000..53b6d7eb4 --- /dev/null +++ b/src/finiteArea/finiteArea/fam/vectorFamDiv.H @@ -0,0 +1,74 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Typedef + Foam::fam + +Description + Specialisation of fam div for a flux. + +SourceFiles + vectorFamDiv.C + +\*---------------------------------------------------------------------------*/ + +#ifndef vectorFamDiv_H +#define vectorFamDiv_H + +#include "famDiv.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace fam +{ + +/*---------------------------------------------------------------------------*\ + Namespace fam functions Declaration +\*---------------------------------------------------------------------------*/ + +template<> +tmp > +div +( + const edgeScalarField& flux, + const GeometricField& vf, + const word& name +); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C b/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C index 6ca7b9db2..c9fef3221 100644 --- a/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C +++ b/src/finiteArea/finiteArea/gradSchemes/gaussFaGrad/gaussFaGrad.C @@ -66,12 +66,7 @@ gaussGrad::grad GeometricField& gGrad = tgGrad(); - gGrad -= vsf*fac::edgeIntegrate(vsf.mesh().Le()); - - // Remove component of gradient normal to surface (area) - const areaVectorField& n = vsf.mesh().faceAreaNormals(); - - gGrad -= n*(n & gGrad); + // Removed for consistencty. Matthias Rauter, 6/Dec/2016 gGrad.correctBoundaryConditions(); gGrad.rename("grad(" + vsf.name() + ')'); @@ -96,15 +91,9 @@ void gaussGrad::correctBoundaryConditions if (!vsf.boundaryField()[patchI].coupled()) { vectorField m = - vsf.mesh().Le().boundaryField()[patchI] - /vsf.mesh().magLe().boundaryField()[patchI]; + vsf.mesh().Le().boundaryField()[patchI]/ + vsf.mesh().magLe().boundaryField()[patchI]; - // Zeljko Tukovic -// gGrad.boundaryField()[patchI] = -// m*vsf.boundaryField()[patchI].snGrad(); - - //HJ Not sure: should this be a complete correction or just the - // tangential part? HJ, 24/Jul/2009 gGrad.boundaryField()[patchI] += m* ( vsf.boundaryField()[patchI].snGrad() diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C index d9796d613..a910159e3 100644 --- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C +++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C @@ -138,11 +138,8 @@ leastSquaresFaGrad::grad } } - // Remove component of gradient normal to surface (area) - const areaVectorField& n = vsf.mesh().faceAreaNormals(); - - lsGrad -= n*(n & lsGrad); + // Removed for consistencty. Matthias Rauter, 6/Dec/2016 lsGrad.correctBoundaryConditions(); gaussGrad::correctBoundaryConditions(vsf, lsGrad); diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C index ed907ca08..09dfb6e25 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.C @@ -62,7 +62,7 @@ void Foam::MRFZone::setMRFFaces() if (cellZoneID_ != -1) { const labelList& cellLabels = mesh_.cellZones()[cellZoneID_]; - forAll(cellLabels, i) + forAll (cellLabels, i) { zoneCell[cellLabels[i]] = true; } @@ -84,15 +84,15 @@ void Foam::MRFZone::setMRFFaces() labelHashSet excludedPatches(excludedPatchLabels_); - forAll(patches, patchI) + forAll (patches, patchI) { const polyPatch& pp = patches[patchI]; if (pp.isWall() && !excludedPatches.found(patchI)) { - forAll(pp, i) + forAll (pp, i) { - label faceI = pp.start()+i; + label faceI = pp.start() + i; if (zoneCell[own[faceI]]) { @@ -103,9 +103,9 @@ void Foam::MRFZone::setMRFFaces() } else if (!isA(pp)) { - forAll(pp, i) + forAll (pp, i) { - label faceI = pp.start()+i; + label faceI = pp.start() + i; if (zoneCell[own[faceI]]) { @@ -138,11 +138,11 @@ void Foam::MRFZone::setMRFFaces() labelList nIncludedFaces(patches.size(), 0); labelList nExcludedFaces(patches.size(), 0); - forAll(patches, patchi) + forAll (patches, patchi) { const polyPatch& pp = patches[patchi]; - forAll(pp, patchFacei) + forAll (pp, patchFacei) { label faceI = pp.start()+patchFacei; @@ -159,7 +159,7 @@ void Foam::MRFZone::setMRFFaces() includedFaces_.setSize(patches.size()); excludedFaces_.setSize(patches.size()); - forAll(nIncludedFaces, patchi) + forAll (nIncludedFaces, patchi) { includedFaces_[patchi].setSize(nIncludedFaces[patchi]); excludedFaces_[patchi].setSize(nExcludedFaces[patchi]); @@ -167,11 +167,11 @@ void Foam::MRFZone::setMRFFaces() nIncludedFaces = 0; nExcludedFaces = 0; - forAll(patches, patchi) + forAll (patches, patchi) { const polyPatch& pp = patches[patchi]; - forAll(pp, patchFacei) + forAll (pp, patchFacei) { label faceI = pp.start()+patchFacei; @@ -201,9 +201,9 @@ void Foam::MRFZone::setMRFFaces() internalFaces.write(); faceSet MRFFaces(mesh_, "includedFaces", 100); - forAll(includedFaces_, patchi) + forAll (includedFaces_, patchi) { - forAll(includedFaces_[patchi], i) + forAll (includedFaces_[patchi], i) { label patchFacei = includedFaces_[patchi][i]; MRFFaces.insert(patches[patchi].start()+patchFacei); @@ -215,12 +215,12 @@ void Foam::MRFZone::setMRFFaces() MRFFaces.write(); faceSet excludedFaces(mesh_, "excludedFaces", 100); - forAll(excludedFaces_, patchi) + forAll (excludedFaces_, patchi) { - forAll(excludedFaces_[patchi], i) + forAll (excludedFaces_[patchi], i) { label patchFacei = excludedFaces_[patchi][i]; - excludedFaces.insert(patches[patchi].start()+patchFacei); + excludedFaces.insert(patches[patchi].start() + patchFacei); } } Pout<< "Writing " << excludedFaces.size() @@ -231,12 +231,30 @@ void Foam::MRFZone::setMRFFaces() } +Foam::vector Foam::MRFZone::Omega() const +{ + if (rampTime_ < SMALL) + { + return omega_.value()*axis_.value(); + } + else + { + // Ramping + Info<< "MRF Ramp: " << Foam::min(mesh_.time().value()/rampTime_, 1) + << endl; + + return Foam::min(mesh_.time().value()/rampTime_, 1)* + omega_.value()*axis_.value(); + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is) : - mesh_(mesh), name_(is), + mesh_(mesh), dict_(is), cellZoneID_(mesh_.cellZones().findZoneID(name_)), excludedPatchNames_ @@ -246,7 +264,7 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is) origin_(dict_.lookup("origin")), axis_(dict_.lookup("axis")), omega_(dict_.lookup("omega")), - Omega_("Omega", omega_*axis_) + rampTime_(dict_.lookupOrDefault("rampTime", 0)) { if (dict_.found("patches")) { @@ -260,12 +278,19 @@ Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is) const polyBoundaryMesh& patches = mesh_.boundaryMesh(); + if (mag(axis_.value()) < SMALL) + { + FatalErrorIn("MRFZone(const fvMesh&, Istream&)") + << "Axis vector has zero magnitude: " << axis_ + << ". This is not allowed" + << abort(FatalError); + } + axis_ = axis_/mag(axis_); - Omega_ = omega_*axis_; excludedPatchLabels_.setSize(excludedPatchNames_.size()); - forAll(excludedPatchNames_, i) + forAll (excludedPatchNames_, i) { excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]); @@ -309,12 +334,12 @@ void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const const scalarField& V = mesh_.V(); vectorField& Usource = UEqn.source(); const vectorField& U = UEqn.psi(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); - forAll(cells, i) + forAll (cells, i) { label celli = cells[i]; - Usource[celli] -= V[celli]*(Omega ^ U[celli]); + Usource[celli] -= V[celli]*(rotVel ^ U[celli]); } } @@ -334,12 +359,12 @@ void Foam::MRFZone::addCoriolis const scalarField& V = mesh_.V(); vectorField& Usource = UEqn.source(); const vectorField& U = UEqn.psi(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); - forAll(cells, i) + forAll (cells, i) { label celli = cells[i]; - Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]); + Usource[celli] -= V[celli]*rho[celli]*(rotVel ^ U[celli]); } } @@ -349,20 +374,20 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const const volVectorField& C = mesh_.C(); const vector& origin = origin_.value(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); const labelList& cells = mesh_.cellZones()[cellZoneID_]; - forAll(cells, i) + forAll (cells, i) { label celli = cells[i]; - U[celli] -= (Omega ^ (C[celli] - origin)); + U[celli] -= (rotVel ^ (C[celli] - origin)); } // Included patches - forAll(includedFaces_, patchi) + forAll (includedFaces_, patchi) { - forAll(includedFaces_[patchi], i) + forAll (includedFaces_[patchi], i) { label patchFacei = includedFaces_[patchi][i]; U.boundaryField()[patchi][patchFacei] = vector::zero; @@ -370,13 +395,13 @@ void Foam::MRFZone::relativeVelocity(volVectorField& U) const } // Excluded patches - forAll(excludedFaces_, patchi) + forAll (excludedFaces_, patchi) { - forAll(excludedFaces_[patchi], i) + forAll (excludedFaces_[patchi], i) { label patchFacei = excludedFaces_[patchi][i]; U.boundaryField()[patchi][patchFacei] -= - (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin)); + (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin)); } } } @@ -387,35 +412,35 @@ void Foam::MRFZone::absoluteVelocity(volVectorField& U) const const volVectorField& C = mesh_.C(); const vector& origin = origin_.value(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); const labelList& cells = mesh_.cellZones()[cellZoneID_]; - forAll(cells, i) + forAll (cells, i) { label celli = cells[i]; - U[celli] += (Omega ^ (C[celli] - origin)); + U[celli] += (rotVel ^ (C[celli] - origin)); } // Included patches - forAll(includedFaces_, patchi) + forAll (includedFaces_, patchi) { - forAll(includedFaces_[patchi], i) + forAll (includedFaces_[patchi], i) { label patchFacei = includedFaces_[patchi][i]; U.boundaryField()[patchi][patchFacei] = - (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin)); + (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin)); } } // Excluded patches - forAll(excludedFaces_, patchi) + forAll (excludedFaces_, patchi) { - forAll(excludedFaces_[patchi], i) + forAll (excludedFaces_[patchi], i) { label patchFacei = excludedFaces_[patchi][i]; U.boundaryField()[patchi][patchFacei] += - (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin)); + (rotVel ^ (C.boundaryField()[patchi][patchFacei] - origin)); } } } @@ -460,36 +485,36 @@ void Foam::MRFZone::meshPhi const surfaceVectorField& Cf = mesh_.Cf(); const vector& origin = origin_.value(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); // Internal faces - forAll(internalFaces_, i) + forAll (internalFaces_, i) { label facei = internalFaces_[i]; - phi[facei] = (Omega ^ (Cf[facei] - origin)); + phi[facei] = (rotVel ^ (Cf[facei] - origin)); } // Included patches - forAll(includedFaces_, patchi) + forAll (includedFaces_, patchi) { - forAll(includedFaces_[patchi], i) + forAll (includedFaces_[patchi], i) { label patchFacei = includedFaces_[patchi][i]; phi.boundaryField()[patchi][patchFacei] = - (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin)); + (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)); } } // Excluded patches - forAll(excludedFaces_, patchi) + forAll (excludedFaces_, patchi) { - forAll(excludedFaces_[patchi], i) + forAll (excludedFaces_[patchi], i) { label patchFacei = excludedFaces_[patchi][i]; phi.boundaryField()[patchi][patchFacei] = - (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin)); + (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)); } } } @@ -497,20 +522,20 @@ void Foam::MRFZone::meshPhi void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const { const vector& origin = origin_.value(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); // Included patches - forAll(includedFaces_, patchi) + forAll (includedFaces_, patchi) { const vectorField& patchC = mesh_.Cf().boundaryField()[patchi]; vectorField pfld(U.boundaryField()[patchi]); - forAll(includedFaces_[patchi], i) + forAll (includedFaces_[patchi], i) { label patchFacei = includedFaces_[patchi][i]; - pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin)); + pfld[patchFacei] = (rotVel ^ (patchC[patchFacei] - origin)); } U.boundaryField()[patchi] == pfld; @@ -529,12 +554,12 @@ void Foam::MRFZone::Su const volVectorField& C = mesh_.C(); const vector& origin = origin_.value(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); - forAll(cells, i) + forAll (cells, i) { source[cells[i]] = - (Omega ^ (C[cells[i]] - origin)) & gradPhi[cells[i]]; + (rotVel ^ (C[cells[i]] - origin)) & gradPhi[cells[i]]; } } @@ -550,13 +575,13 @@ void Foam::MRFZone::Su const volVectorField& C = mesh_.C(); const vector& origin = origin_.value(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); - forAll(cells, i) + forAll (cells, i) { source[cells[i]] = - ((Omega ^ (C[cells[i]] - origin)) & gradPhi[cells[i]]) - - (Omega ^ phi[cells[i]]); + ((rotVel ^ (C[cells[i]] - origin)) & gradPhi[cells[i]]) + - (rotVel ^ phi[cells[i]]); } } diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H index 63cd0dd16..954bc85b1 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZone.H +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZone.H @@ -29,7 +29,7 @@ Description obtained from a control dictionary constructed from the given stream. The rotation of the MRF region is defined by an origin and axis of - rotation and an angular speed. + rotation and an angular speed in rad/s. SourceFiles MRFZone.C @@ -65,15 +65,22 @@ class MRFZone { // Private data - const fvMesh& mesh_; - + //- Name of the MRF zone const word name_; + //- Mesh reference + const fvMesh& mesh_; + + //- Reference to dictionary const dictionary dict_; + //- MRF cell zone ID: cells in MRF region label cellZoneID_; + //- List of patch names touching the MRF zone which do not move const wordList excludedPatchNames_; + + //- List of patch labels touching the MRF zone which do not move labelList excludedPatchLabels_; //- Internal faces that are part of MRF @@ -85,17 +92,34 @@ class MRFZone //- Excluded faces (per patch) that do not move with the MRF labelListList excludedFaces_; + //- Origin of rotation const dimensionedVector origin_; + + //- Axis of rotation dimensionedVector axis_; - const dimensionedScalar omega_; - dimensionedVector Omega_; + + //- Rotational velocity in rad/s + dimensionedScalar omega_; + + //- Ramping time scale + scalar rampTime_; // Private Member Functions + //- Disallow default bitwise copy construct + MRFZone(const MRFZone&); + + //- Disallow default bitwise assignment + void operator=(const MRFZone&); + + //- Divide faces in frame according to patch void setMRFFaces(); + //- Return rotational vector + vector Omega() const; + //- Make the given absolute mass/vol flux relative within MRF region template void relativeRhoFlux @@ -112,12 +136,6 @@ class MRFZone surfaceScalarField& phi ) const; - //- Disallow default bitwise copy construct - MRFZone(const MRFZone&); - - //- Disallow default bitwise assignment - void operator=(const MRFZone&); - public: @@ -170,7 +188,11 @@ public: void addCoriolis(fvVectorMatrix& UEqn) const; //- Add the Coriolis force contribution to the momentum equation - void addCoriolis(const volScalarField& rho, fvVectorMatrix& UEqn) const; + void addCoriolis + ( + const volScalarField& rho, + fvVectorMatrix& UEqn + ) const; //- Make the given absolute velocity relative within the MRF region void relativeVelocity(volVectorField& U) const; diff --git a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C index d8070148b..0ae07364a 100644 --- a/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C +++ b/src/finiteVolume/cfdTools/general/MRF/MRFZoneTemplates.C @@ -41,19 +41,19 @@ void Foam::MRFZone::relativeRhoFlux const surfaceVectorField& Sf = mesh_.Sf(); const vector& origin = origin_.value(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); // Internal faces - forAll(internalFaces_, i) + forAll (internalFaces_, i) { label facei = internalFaces_[i]; - phi[facei] -= rho[facei]*(Omega ^ (Cf[facei] - origin)) & Sf[facei]; + phi[facei] -= rho[facei]*(rotVel ^ (Cf[facei] - origin)) & Sf[facei]; } // Included patches - forAll(includedFaces_, patchi) + forAll (includedFaces_, patchi) { - forAll(includedFaces_[patchi], i) + forAll (includedFaces_[patchi], i) { label patchFacei = includedFaces_[patchi][i]; @@ -62,15 +62,15 @@ void Foam::MRFZone::relativeRhoFlux } // Excluded patches - forAll(excludedFaces_, patchi) + forAll (excludedFaces_, patchi) { - forAll(excludedFaces_[patchi], i) + forAll (excludedFaces_[patchi], i) { label patchFacei = excludedFaces_[patchi][i]; phi.boundaryField()[patchi][patchFacei] -= rho.boundaryField()[patchi][patchFacei] - *(Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) + *(rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) & Sf.boundaryField()[patchi][patchFacei]; } } @@ -88,37 +88,37 @@ void Foam::MRFZone::absoluteRhoFlux const surfaceVectorField& Sf = mesh_.Sf(); const vector& origin = origin_.value(); - const vector& Omega = Omega_.value(); + const vector rotVel = Omega(); // Internal faces - forAll(internalFaces_, i) + forAll (internalFaces_, i) { label facei = internalFaces_[i]; - phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei]; + phi[facei] += (rotVel ^ (Cf[facei] - origin)) & Sf[facei]; } // Included patches - forAll(includedFaces_, patchi) + forAll (includedFaces_, patchi) { - forAll(includedFaces_[patchi], i) + forAll (includedFaces_[patchi], i) { label patchFacei = includedFaces_[patchi][i]; phi.boundaryField()[patchi][patchFacei] += - (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) + (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) & Sf.boundaryField()[patchi][patchFacei]; } } // Excluded patches - forAll(excludedFaces_, patchi) + forAll (excludedFaces_, patchi) { - forAll(excludedFaces_[patchi], i) + forAll (excludedFaces_[patchi], i) { label patchFacei = excludedFaces_[patchi][i]; phi.boundaryField()[patchi][patchFacei] += - (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) + (rotVel ^ (Cf.boundaryField()[patchi][patchFacei] - origin)) & Sf.boundaryField()[patchi][patchFacei]; } } diff --git a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C index 25f33159b..a0dca512b 100644 --- a/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C +++ b/src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C @@ -62,7 +62,7 @@ void Foam::setRefCell ")", dict ) << "Illegal master cellID " << refCelli - << ". Should be 0.." << field.mesh().nCells() + << ". Should be 0." << field.mesh().nCells() << exit(FatalIOError); } } diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C index 0a92f62b1..4d6d58a1d 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.C @@ -31,27 +31,28 @@ License #include "coeffFields.H" #include "demandDrivenData.H" -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template -processorFvPatchField::processorFvPatchField +Foam::processorFvPatchField::processorFvPatchField ( const fvPatch& p, const DimensionedField& iF ) : coupledFvPatchField(p, iF), - procPatch_(refCast(p)) + procPatch_(refCast(p)), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + sendBuf_(0), + receiveBuf_(0), + scalarSendBuf_(0), + scalarReceiveBuf_(0) {} template -processorFvPatchField::processorFvPatchField +Foam::processorFvPatchField::processorFvPatchField ( const fvPatch& p, const DimensionedField& iF, @@ -59,13 +60,18 @@ processorFvPatchField::processorFvPatchField ) : coupledFvPatchField(p, iF, f), - procPatch_(refCast(p)) + procPatch_(refCast(p)), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + sendBuf_(0), + receiveBuf_(0), + scalarSendBuf_(0), + scalarReceiveBuf_(0) {} -// Construct by mapping given processorFvPatchField template -processorFvPatchField::processorFvPatchField +Foam::processorFvPatchField::processorFvPatchField ( const processorFvPatchField& ptf, const fvPatch& p, @@ -74,9 +80,15 @@ processorFvPatchField::processorFvPatchField ) : coupledFvPatchField(ptf, p, iF, mapper), - procPatch_(refCast(p)) + procPatch_(refCast(p)), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + sendBuf_(0), + receiveBuf_(0), + scalarSendBuf_(0), + scalarReceiveBuf_(0) { - if (!isType(this->patch())) + if (!isA(this->patch())) { FatalErrorIn ( @@ -98,7 +110,7 @@ processorFvPatchField::processorFvPatchField template -processorFvPatchField::processorFvPatchField +Foam::processorFvPatchField::processorFvPatchField ( const fvPatch& p, const DimensionedField& iF, @@ -106,9 +118,15 @@ processorFvPatchField::processorFvPatchField ) : coupledFvPatchField(p, iF, dict), - procPatch_(refCast(p)) + procPatch_(refCast(p)), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + sendBuf_(0), + receiveBuf_(0), + scalarSendBuf_(0), + scalarReceiveBuf_(0) { - if (!isType(p)) + if (!isA(p)) { FatalIOErrorIn ( @@ -130,42 +148,79 @@ processorFvPatchField::processorFvPatchField template -processorFvPatchField::processorFvPatchField +Foam::processorFvPatchField::processorFvPatchField ( const processorFvPatchField& ptf ) : processorLduInterfaceField(), coupledFvPatchField(ptf), - procPatch_(refCast(ptf.patch())) + procPatch_(refCast(ptf.patch())), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + sendBuf_(0), + receiveBuf_(0), + scalarSendBuf_(0), + scalarReceiveBuf_(0) {} template -processorFvPatchField::processorFvPatchField +Foam::processorFvPatchField::processorFvPatchField ( const processorFvPatchField& ptf, const DimensionedField& iF ) : coupledFvPatchField(ptf, iF), - procPatch_(refCast(ptf.patch())) -{} + procPatch_(refCast(ptf.patch())), + outstandingSendRequest_(-1), + outstandingRecvRequest_(-1), + sendBuf_(0), + receiveBuf_(0), + scalarSendBuf_(0), + scalarReceiveBuf_(0) +{ + if (debug && !ptf.ready()) + { + FatalErrorIn + ( + "processorFvPatchField::processorFvPatchField\n" + "(\n" + " const processorFvPatchField& ptf,\n" + " const DimensionedField& iF\n" + ")" + ) << "On patch " << procPatch_.name() << " outstanding request." + << abort(FatalError); + } +} // * * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * // template -processorFvPatchField::~processorFvPatchField() +Foam::processorFvPatchField::~processorFvPatchField() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template -tmp > processorFvPatchField::patchNeighbourField() const +Foam::tmp > +Foam::processorFvPatchField::patchNeighbourField() const { - // Warning: returning own patch field, which on update stores + if (debug && !this->ready()) + { + FatalErrorIn + ( + "tmp >" + "processorFvPatchField::patchNeighbourField() const" + ) << "On patch " << procPatch_.name() + << " outstanding request." + << abort(FatalError); + } + + // Warning: returning own patch field, which only after update stores // actual neighbour data // HJ, 14/May/2009 return *this; @@ -173,27 +228,79 @@ tmp > processorFvPatchField::patchNeighbourField() const template -void processorFvPatchField::initEvaluate +void Foam::processorFvPatchField::initEvaluate ( const Pstream::commsTypes commsType ) { if (Pstream::parRun()) { - procPatch_.compressedSend(commsType, this->patchInternalField()()); + // Collect data into send buffer + sendBuf_ = this->patchInternalField(); + + if (commsType == Pstream::nonBlocking) + { + // Fast path. Receive into *this + + outstandingRecvRequest_ = Pstream::nRequests(); + + IPstream::read + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(this->begin()), + this->byteSize(), + procPatch_.tag(), + procPatch_.comm() + ); + + outstandingSendRequest_ = Pstream::nRequests(); + + OPstream::write + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(sendBuf_.begin()), + this->byteSize(), + procPatch_.tag(), + procPatch_.comm() + ); + } + else + { + procPatch_.send(commsType, sendBuf_); + } } } template -void processorFvPatchField::evaluate +void Foam::processorFvPatchField::evaluate ( const Pstream::commsTypes commsType ) { if (Pstream::parRun()) { - procPatch_.compressedReceive(commsType, *this); + if (commsType == Pstream::nonBlocking) + { + // Fast path. Received into *this + + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + Pstream::waitRequest(outstandingRecvRequest_); + } + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + } + else + { + procPatch_.receive(commsType, *this); + } if (doTransform()) { @@ -204,14 +311,15 @@ void processorFvPatchField::evaluate template -tmp > processorFvPatchField::snGrad() const +Foam::tmp > +Foam::processorFvPatchField::snGrad() const { return this->patch().deltaCoeffs()*(*this - this->patchInternalField()); } template -void processorFvPatchField::initInterfaceMatrixUpdate +void Foam::processorFvPatchField::initInterfaceMatrixUpdate ( const scalarField& psiInternal, scalarField&, @@ -222,16 +330,68 @@ void processorFvPatchField::initInterfaceMatrixUpdate const bool switchToLhs ) const { - procPatch_.compressedSend - ( - commsType, - this->patch().patchInternalField(psiInternal)() - ); + scalarSendBuf_ = this->patch().patchInternalField(psiInternal); + + if (commsType == Pstream::nonBlocking) + { + // Fast path. + if (debug && !this->ready()) + { + FatalErrorIn + ( + "void processorFvPatchField::initInterfaceMatrixUpdate\n" + "(\n" + " const scalarField& psiInternal,\n" + " scalarField&,\n" + " const lduMatrix&,\n" + " const scalarField&,\n" + " const direction,\n" + " const Pstream::commsTypes commsType,\n" + " const bool switchToLhs\n" + ") const" + ) << "On patch " << procPatch_.name() + << " outstanding request." + << abort(FatalError); + } + + // Check buffer size + scalarReceiveBuf_.setSize(this->size()); + + outstandingRecvRequest_ = Pstream::nRequests(); + + IPstream::read + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(scalarReceiveBuf_.begin()), + scalarReceiveBuf_.byteSize(), + procPatch_.tag(), + procPatch_.comm() + ); + + outstandingSendRequest_ = Pstream::nRequests(); + + OPstream::write + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(scalarSendBuf_.begin()), + scalarSendBuf_.byteSize(), + procPatch_.tag(), + procPatch_.comm() + ); + } + else + { + procPatch_.send(commsType, scalarSendBuf_); + } + + const_cast&>(*this).updatedMatrix() = false; } template -void processorFvPatchField::updateInterfaceMatrix +void Foam::processorFvPatchField::updateInterfaceMatrix ( const scalarField&, scalarField& result, @@ -242,36 +402,65 @@ void processorFvPatchField::updateInterfaceMatrix const bool switchToLhs ) const { - scalarField pnf - ( - procPatch_.compressedReceive(commsType, this->size())() - ); + if (this->updatedMatrix()) + { + return; + } + + if (commsType == Pstream::nonBlocking) + { + // Fast path. + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + Pstream::waitRequest(outstandingRecvRequest_); + } + + // Recv finished so assume sending finished as well. + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + } + else + { + // Check size + scalarReceiveBuf_.setSize(this->size()); + + procPatch_.receive(commsType, scalarReceiveBuf_); + } + + // The data is now in scalarReceiveBuf_ for both cases // Transform according to the transformation tensor - transformCoupleField(pnf, cmpt); + transformCoupleField(scalarReceiveBuf_, cmpt); // Multiply the field by coefficients and add into the result + const unallocLabelList& faceCells = this->patch().faceCells(); if (switchToLhs) { forAll(faceCells, elemI) { - result[faceCells[elemI]] += coeffs[elemI]*pnf[elemI]; + result[faceCells[elemI]] += coeffs[elemI]*scalarReceiveBuf_[elemI]; } } else { forAll(faceCells, elemI) { - result[faceCells[elemI]] -= coeffs[elemI]*pnf[elemI]; + result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI]; } } + + const_cast&>(*this).updatedMatrix() = true; } template -void processorFvPatchField::initInterfaceMatrixUpdate +void Foam::processorFvPatchField::initInterfaceMatrixUpdate ( const Field& psiInternal, Field&, @@ -281,16 +470,71 @@ void processorFvPatchField::initInterfaceMatrixUpdate const bool switchToLhs ) const { - procPatch_.compressedSend - ( - commsType, - this->patch().patchInternalField(psiInternal)() - ); + sendBuf_ = this->patch().patchInternalField(psiInternal); + + if (commsType == Pstream::nonBlocking) + { + // Fast path. + if (debug && !this->ready()) + { + FatalErrorIn + ( + "void processorFvPatchField::initInterfaceMatrixUpdate\n" + "(\n" + " const Field& psiInternal,\n" + " Field&,\n" + " const BlockLduMatrix&,\n" + " const CoeffField&,\n" + " const Pstream::commsTypes commsType,\n" + " const bool switchToLhs\n" + ") const" + ) << "On patch " << procPatch_.name() + << " outstanding request." + << abort(FatalError); + } + + // Check buffer size + receiveBuf_.setSize(this->size()); + + outstandingRecvRequest_ = Pstream::nRequests(); + + IPstream::read + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(receiveBuf_.begin()), + receiveBuf_.byteSize(), + procPatch_.tag(), + procPatch_.comm() + ); + + outstandingSendRequest_ = Pstream::nRequests(); + + OPstream::write + ( + Pstream::nonBlocking, + procPatch_.neighbProcNo(), + reinterpret_cast(sendBuf_.begin()), + sendBuf_.byteSize(), + procPatch_.tag(), + procPatch_.comm() + ); + } + else + { + procPatch_.send + ( + commsType, + this->patch().patchInternalField(psiInternal)() + ); + } + + const_cast&>(*this).updatedMatrix() = false; } template -void processorFvPatchField::updateInterfaceMatrix +void Foam::processorFvPatchField::updateInterfaceMatrix ( const Field& psiInternal, Field& result, @@ -300,14 +544,40 @@ void processorFvPatchField::updateInterfaceMatrix const bool switchToLhs ) const { - Field pnf - ( - procPatch_.compressedReceive(commsType, this->size())() - ); + if (this->updatedMatrix()) + { + return; + } - // Multiply neighbour field with coeffs and re-use pnf for result + if (commsType == Pstream::nonBlocking) + { + // Fast path. + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + Pstream::waitRequest(outstandingRecvRequest_); + } + + // Recv finished so assume sending finished as well. + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + } + else + { + // Check size + receiveBuf_.setSize(this->size()); + + procPatch_.receive(commsType, receiveBuf_); + } + + // The data is now in receiveBuf_ for both cases + + // Multiply neighbour field with coeffs and re-use buffer for result // of multiplication - multiply(pnf, coeffs, pnf); + multiply(receiveBuf_, coeffs, receiveBuf_); const unallocLabelList& faceCells = this->patch().faceCells(); @@ -315,21 +585,56 @@ void processorFvPatchField::updateInterfaceMatrix { forAll(faceCells, elemI) { - result[faceCells[elemI]] += pnf[elemI]; + result[faceCells[elemI]] += receiveBuf_[elemI]; } } else { forAll(faceCells, elemI) { - result[faceCells[elemI]] -= pnf[elemI]; + result[faceCells[elemI]] -= receiveBuf_[elemI]; } } + + const_cast&>(*this).updatedMatrix() = true; } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +template +bool Foam::processorFvPatchField::ready() const +{ + if + ( + outstandingSendRequest_ >= 0 + && outstandingSendRequest_ < Pstream::nRequests() + ) + { + bool finished = Pstream::finishedRequest(outstandingSendRequest_); + + if (!finished) + { + return false; + } + } + outstandingSendRequest_ = -1; + + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + bool finished = Pstream::finishedRequest(outstandingRecvRequest_); + + if (!finished) + { + return false; + } + } + outstandingRecvRequest_ = -1; + + return true; +} -} // End namespace Foam // ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H index 4f23bcc25..05972d2ae 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchField.H @@ -69,6 +69,25 @@ class processorFvPatchField //- Local reference cast into the processor patch const processorFvPatch& procPatch_; + // Sending and receiving + + //- Outstanding request + mutable label outstandingSendRequest_; + + //- Outstanding request + mutable label outstandingRecvRequest_; + + //- Send buffer. + mutable Field sendBuf_; + + //- Receive buffer. + mutable Field receiveBuf_; + + //- Scalar send buffer + mutable Field scalarSendBuf_; + + //- Scalar receive buffer + mutable Field scalarReceiveBuf_; public: @@ -143,7 +162,7 @@ public: //- Destructor - ~processorFvPatchField(); + virtual ~processorFvPatchField(); // Member functions @@ -249,7 +268,7 @@ public: ) const; - //- Processor coupled interface functions + // Processor coupled interface functions //- Return processor number virtual int myProcNo() const @@ -263,6 +282,19 @@ public: return procPatch_.neighbProcNo(); } + + // Communication support + + //- Is all data available + virtual bool ready() const; + + //- Return communicator used for parallel communication + virtual int comm() const + { + return procPatch_.comm(); + } + + //- Does the patch field perform the transfromation virtual bool doTransform() const { diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C index 105c6b466..acce48d51 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.C @@ -32,26 +32,6 @@ namespace Foam // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template<> -void processorFvPatchField::initInterfaceMatrixUpdate -( - const scalarField& psiInternal, - scalarField&, - const lduMatrix&, - const scalarField&, - const direction, - const Pstream::commsTypes commsType, - const bool switchToLhs -) const -{ - procPatch_.compressedSend - ( - commsType, - patch().patchInternalField(psiInternal)() - ); -} - - template<> void processorFvPatchField::updateInterfaceMatrix ( @@ -64,82 +44,60 @@ void processorFvPatchField::updateInterfaceMatrix const bool switchToLhs ) const { - scalarField pnf - ( - procPatch_.compressedReceive(commsType, this->size())() - ); + if (this->updatedMatrix()) + { + return; + } - const unallocLabelList& faceCells = patch().faceCells(); + if (commsType == Pstream::nonBlocking) + { + // Fast path. + if + ( + outstandingRecvRequest_ >= 0 + && outstandingRecvRequest_ < Pstream::nRequests() + ) + { + Pstream::waitRequest(outstandingRecvRequest_); + } + + // Recv finished so assume sending finished as well. + outstandingSendRequest_ = -1; + outstandingRecvRequest_ = -1; + } + else + { + // Check size + scalarReceiveBuf_.setSize(this->size()); + + procPatch_.receive(commsType, scalarReceiveBuf_); + } + + // The data is now in scalarReceiveBuf_ for both cases + + // Transform according to the transformation tensor + // No transform for scalar. HJ, 29/Nov/2016 + + // Multiply the field by coefficients and add into the result + + const unallocLabelList& faceCells = this->patch().faceCells(); if (switchToLhs) { - forAll (faceCells, facei) + forAll(faceCells, elemI) { - result[faceCells[facei]] += coeffs[facei]*pnf[facei]; + result[faceCells[elemI]] += coeffs[elemI]*scalarReceiveBuf_[elemI]; } } else { - forAll (faceCells, facei) + forAll(faceCells, elemI) { - result[faceCells[facei]] -= coeffs[facei]*pnf[facei]; + result[faceCells[elemI]] -= coeffs[elemI]*scalarReceiveBuf_[elemI]; } } -} - -template<> -void processorFvPatchField::initInterfaceMatrixUpdate -( - const scalarField& psiInternal, - scalarField&, - const BlockLduMatrix&, - const CoeffField&, - const Pstream::commsTypes commsType, - const bool switchToLhs -) const -{ - procPatch_.compressedSend - ( - commsType, - patch().patchInternalField(psiInternal)() - ); -} - - -template<> -void processorFvPatchField::updateInterfaceMatrix -( - const scalarField&, - scalarField& result, - const BlockLduMatrix&, - const CoeffField& coeffs, - const Pstream::commsTypes commsType, - const bool switchToLhs -) const -{ - scalarField pnf - ( - procPatch_.compressedReceive(commsType, this->size())() - ); - - const unallocLabelList& faceCells = patch().faceCells(); - const scalarField& scalarCoeffs = coeffs.asScalar(); - - if (switchToLhs) - { - forAll (faceCells, facei) - { - result[faceCells[facei]] += scalarCoeffs[facei]*pnf[facei]; - } - } - else - { - forAll (faceCells, facei) - { - result[faceCells[facei]] -= scalarCoeffs[facei]*pnf[facei]; - } - } + const_cast&>(*this).updatedMatrix() = true; } diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.H b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.H index 4c8e21a54..1c9d84b79 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/processor/processorFvPatchScalarField.H @@ -36,18 +36,8 @@ namespace Foam // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -template<> -void processorFvPatchField::initInterfaceMatrixUpdate -( - const scalarField&, - scalarField&, - const lduMatrix&, - const scalarField&, - const direction, - const Pstream::commsTypes commsType, - const bool switchToLhs -) const; - +// Removed unnecessary specialisation for initInterfaceMatrixUpdate +// HJ, 29/Nov/2016 template<> void processorFvPatchField::updateInterfaceMatrix @@ -62,29 +52,11 @@ void processorFvPatchField::updateInterfaceMatrix ) const; -template<> -void processorFvPatchField::initInterfaceMatrixUpdate -( - const scalarField& psiInternal, - scalarField&, - const BlockLduMatrix&, - const CoeffField&, - const Pstream::commsTypes commsType, - const bool switchToLhs -) const; +// Removed unnecessary specialisation for initInterfaceMatrixUpdate +// and updateInterfaceMatrix for block matrices +// HJ, 29/Nov/2016 -template<> -void processorFvPatchField::updateInterfaceMatrix -( - const scalarField&, - scalarField& result, - const BlockLduMatrix&, - const CoeffField& coeffs, - const Pstream::commsTypes commsType, - const bool switchToLhs -) const; - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.C index f714fc01e..54b6d50e9 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.C @@ -141,7 +141,10 @@ tmp > wedgeFvPatchField::snGrad() const Field pif = this->patchInternalField(); return ( - transform(refCast(this->patch()).cellT(), pif) - pif + transform + ( + refCast(this->patch()).cellT(), pif + ) - pif )*(0.5*this->patch().deltaCoeffs()); } diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.H index 3dce3802d..a4ac779c5 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/wedge/wedgeFvPatchField.H @@ -44,7 +44,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class wedgeFvPatch Declaration + Class wedgeFvPatchField Declaration \*---------------------------------------------------------------------------*/ template diff --git a/src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletOutletVelocity/pressureDirectedInletOutletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletOutletVelocity/pressureDirectedInletOutletVelocityFvPatchVectorField.C index 6d9425903..83b9e1c8b 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletOutletVelocity/pressureDirectedInletOutletVelocityFvPatchVectorField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletOutletVelocity/pressureDirectedInletOutletVelocityFvPatchVectorField.C @@ -46,7 +46,7 @@ pressureDirectedInletOutletVelocityFvPatchVectorField mixedFvPatchVectorField(p, iF), phiName_("phi"), rhoName_("rho"), - inletDir_(p.size()) + inletDir_(p.size(), vector(1, 0, 0)) { refValue() = *this; refGrad() = vector::zero; @@ -84,6 +84,30 @@ pressureDirectedInletOutletVelocityFvPatchVectorField inletDir_("inletDirection", dict, p.size()) { fvPatchVectorField::operator=(vectorField("value", dict, p.size())); + + if (!inletDir_.empty()) + { + if (min(mag(inletDir_)) < SMALL) + { + FatalErrorIn + ( + "pressureDirectedInletOutletVelocityFvPatchVectorField::\n" + "pressureDirectedInletOutletVelocityFvPatchVectorField\n" + "(\n" + " const fvPatch& p,\n" + " const DimensionedField& iF,\n" + " const dictionary& dict\n" + ")" + ) << "Badly defined inlet direction for field " + << this->dimensionedInternalField().name() + << " and patch " << this->patch().name() + << abort(FatalError); + } + } + + // Normalise to obtain the flow direction + inletDir_ /= (mag(inletDir_) + SMALL); + refValue() = *this; refGrad() = vector::zero; valueFraction() = 0.0; diff --git a/src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletVelocity/pressureDirectedInletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletVelocity/pressureDirectedInletVelocityFvPatchVectorField.C index 5b7760545..61b713b96 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletVelocity/pressureDirectedInletVelocityFvPatchVectorField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/pressureDirectedInletVelocity/pressureDirectedInletVelocityFvPatchVectorField.C @@ -46,7 +46,7 @@ pressureDirectedInletVelocityFvPatchVectorField fixedValueFvPatchVectorField(p, iF), phiName_("phi"), rhoName_("rho"), - inletDir_(p.size()) + inletDir_(p.size(), vector(1, 0, 0)) {} @@ -80,6 +80,29 @@ pressureDirectedInletVelocityFvPatchVectorField inletDir_("inletDirection", dict, p.size()) { fvPatchVectorField::operator=(vectorField("value", dict, p.size())); + + if (!inletDir_.empty()) + { + if (min(mag(inletDir_)) < SMALL) + { + FatalErrorIn + ( + "pressureDirectedInletVelocityFvPatchVectorField::\n" + "pressureDirectedInletVelocityFvPatchVectorField\n" + "(\n" + " const fvPatch& p,\n" + " const DimensionedField& iF,\n" + " const dictionary& dict\n" + ")" + ) << "Badly defined inlet direction for field " + << this->dimensionedInternalField().name() + << " and patch " << this->patch().name() + << abort(FatalError); + } + } + + // Normalise to obtain the flow direction + inletDir_ /= (mag(inletDir_) + SMALL); } diff --git a/src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.C index 948de74ac..fddc43131 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/pressureInletOutletVelocity/pressureInletOutletVelocityFvPatchVectorField.C @@ -192,7 +192,6 @@ void pressureInletOutletVelocityFvPatchVectorField::updateCoeffs() } directionMixedFvPatchVectorField::updateCoeffs(); - directionMixedFvPatchVectorField::evaluate(); } diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C index a17e1ff44..793a55d90 100644 --- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C +++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C @@ -698,7 +698,10 @@ void Foam::fvMatrix::relax() { if (psi_.mesh().solutionDict().relaxEquation(psi_.name())) { - relax(psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name())); + relax + ( + psi_.mesh().solutionDict().equationRelaxationFactor(psi_.name()) + ); } else { diff --git a/src/finiteVolume/fvMesh/fvMesh.C b/src/finiteVolume/fvMesh/fvMesh.C index bd69d37e7..928b3b83b 100644 --- a/src/finiteVolume/fvMesh/fvMesh.C +++ b/src/finiteVolume/fvMesh/fvMesh.C @@ -556,6 +556,12 @@ void Foam::fvMesh::updateMesh(const mapPolyMesh& mpm) // This is a temporary solution surfaceInterpolation::movePoints(); + // Note: deltaCoeffs cannot be left on lazy evaluation on mesh motion + // because tangled comms will occur when they are accessed from + // individual boundary conditions + // HJ, VV and IG, 25/Oct/2016 + deltaCoeffs(); + // Function object update moved to polyMesh // HJ, 29/Aug/2010 } @@ -577,6 +583,12 @@ void Foam::fvMesh::syncUpdateMesh() // This is a temporary solution surfaceInterpolation::movePoints(); + // Note: deltaCoeffs cannot be left on lazy evaluation on mesh motion + // because tangled comms will occur when they are accessed from + // individual boundary conditions + // HJ, VV and IG, 25/Oct/2016 + deltaCoeffs(); + // Function object update moved to polyMesh // HJ, 29/Aug/2010 } @@ -681,6 +693,12 @@ Foam::tmp Foam::fvMesh::movePoints(const pointField& p) // Function object update moved to polyMesh // HJ, 29/Aug/2010 + // Note: deltaCoeffs cannot be left on lazy evaluation on mesh motion + // because tangled comms will occur when they are accessed from + // individual boundary conditions + // HJ, VV and IG, 25/Oct/2016 + deltaCoeffs(); + return tsweptVols; } diff --git a/src/finiteVolume/fvMesh/fvMesh.H b/src/finiteVolume/fvMesh/fvMesh.H index 68e9bef2d..86d083a80 100644 --- a/src/finiteVolume/fvMesh/fvMesh.H +++ b/src/finiteVolume/fvMesh/fvMesh.H @@ -123,6 +123,7 @@ class fvMesh // Make geometric data + //- Make face areas void makeSf() const; //- Make face area magnitudes @@ -216,9 +217,8 @@ public: ); - // Destructor - - virtual ~fvMesh(); + //- Destructor + virtual ~fvMesh(); // Member Functions @@ -285,6 +285,15 @@ public: } + // Communication support + + //- Return communicator used for parallel communication + label comm() const + { + return polyMesh::comm(); + } + + //- Return cell volumes const DimensionedField& V() const; diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.H index d47c4079d..35f550a39 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/cyclicGgi/cyclicGgiFvPatch.H @@ -90,10 +90,9 @@ public: {} - // Destructor - - virtual ~cyclicGgiFvPatch() - {} + //- Destructor + virtual ~cyclicGgiFvPatch() + {} // Member functions diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C index 05ad430ce..3b1addbde 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.C @@ -258,6 +258,12 @@ bool Foam::ggiFvPatch::localParallel() const } +const Foam::mapDistribute& Foam::ggiFvPatch::map() const +{ + return ggiPolyPatch_.map(); +} + + const Foam::scalarListList& Foam::ggiFvPatch::weights() const { if (ggiPolyPatch_.master()) @@ -271,6 +277,12 @@ const Foam::scalarListList& Foam::ggiFvPatch::weights() const } +void Foam::ggiFvPatch::expandAddrToZone(labelField& lf) const +{ + lf = ggiPolyPatch_.fastExpand(lf); +} + + Foam::tmp Foam::ggiFvPatch::interfaceInternalField ( const unallocLabelList& internalData @@ -306,6 +318,7 @@ void Foam::ggiFvPatch::initInternalFieldTransfer const unallocLabelList& iF ) const { + // Label transfer is local without global reduction labelTransferBuffer_ = patchInternalField(iF); } diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.H index 6aa101dd3..c063d185d 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/ggi/ggiFvPatch.H @@ -92,13 +92,27 @@ public: {} - // Destructor - - virtual ~ggiFvPatch(); + //- Destructor + virtual ~ggiFvPatch(); // Member functions + // Communication support + + //- Return communicator used for parallel communication + virtual int comm() const + { + return ggiPolyPatch_.comm(); + } + + //- Return message tag used for sending + virtual int tag() const + { + return ggiPolyPatch_.tag(); + } + + // Access //- Return shadow patch @@ -171,6 +185,9 @@ public: //- Is the patch localised on a single processor virtual bool localParallel() const; + //- Return mapDistribute + virtual const mapDistribute& map() const; + //- Return weights. Master side returns own weights and // slave side returns weights from master virtual const scalarListList& weights() const; @@ -187,6 +204,11 @@ public: return coupledFvPatch::reverseT(); } + //- Expand addressing to zone + // Used in optimised AMG coarsening + virtual void expandAddrToZone(labelField&) const; + + //- Return the values of the given internal data adjacent to // the interface as a field virtual tmp interfaceInternalField diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/mixingPlane/mixingPlaneFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/mixingPlane/mixingPlaneFvPatch.H index c03fcab6f..5f4348b4a 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/mixingPlane/mixingPlaneFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/mixingPlane/mixingPlaneFvPatch.H @@ -95,13 +95,27 @@ public: {} - // Destructor - - virtual ~mixingPlaneFvPatch(); + //- Destructor + virtual ~mixingPlaneFvPatch(); // Member functions + // Communication support + + //- Return communicator used for parallel communication + virtual int comm() const + { + return mixingPlanePolyPatch_.comm(); + } + + //- Return message tag used for sending + virtual int tag() const + { + return mixingPlanePolyPatch_.tag(); + } + + // Access //- Return shadow patch diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H index 4e4cbfa6e..77fc4809b 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/processor/processorFvPatch.H @@ -86,10 +86,9 @@ public: {} - // Destructor - - virtual ~processorFvPatch() - {} + //- Destructor + virtual ~processorFvPatch() + {} // Member functions @@ -119,6 +118,21 @@ public: } } + // Communication support + + //- Return communicator used for parallel communication + virtual int comm() const + { + return procPolyPatch_.comm(); + } + + //- Return message tag used for sending + virtual int tag() const + { + return procPolyPatch_.tag(); + } + + //- Return face transformation tensor virtual const tensorField& forwardT() const { diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C index 41e8a58c1..6da84ef02 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.C @@ -34,6 +34,7 @@ Author #include "fvMesh.H" #include "fvBoundaryMesh.H" #include "foamTime.H" +#include "mapDistribute.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -302,6 +303,12 @@ bool Foam::regionCoupleFvPatch::localParallel() const } +const Foam::mapDistribute& Foam::regionCoupleFvPatch::map() const +{ + return rcPolyPatch_.map(); +} + + const Foam::scalarListList& Foam::regionCoupleFvPatch::weights() const { if (rcPolyPatch_.master()) @@ -315,6 +322,16 @@ const Foam::scalarListList& Foam::regionCoupleFvPatch::weights() const } +void Foam::regionCoupleFvPatch::expandAddrToZone(labelField& lf) const +{ + // Missing code. Activate for AMG solvers across regionCoupleFvPatch + notImplemented + ( + "void regionCoupleFvPatch::expandAddrToZone(labelField& lf) const" + ); +} + + Foam::tmp Foam::regionCoupleFvPatch::interfaceInternalField ( const unallocLabelList& internalData diff --git a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H index f09721e88..50b529962 100644 --- a/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H +++ b/src/finiteVolume/fvMesh/fvPatches/constraint/regionCouple/regionCoupleFvPatch.H @@ -95,9 +95,8 @@ public: {} - // Destructor - - virtual ~regionCoupleFvPatch(); + //- Destructor + virtual ~regionCoupleFvPatch(); // Member functions @@ -123,6 +122,21 @@ public: virtual tmp delta() const; + // Communication support + + //- Return communicator used for parallel communication + virtual int comm() const + { + return rcPolyPatch_.comm(); + } + + //- Return message tag used for sending + virtual int tag() const + { + return rcPolyPatch_.tag(); + } + + // Interpolation //- Interpolate face field @@ -180,6 +194,9 @@ public: //- Is the patch localised on a single processor virtual bool localParallel() const; + //- Return mapDistribute + const mapDistribute& map() const; + //- Return weights. Master side returns own weights and // slave side returns weights from master virtual const scalarListList& weights() const; @@ -196,6 +213,11 @@ public: return coupledFvPatch::reverseT(); } + //- Expand addressing to zone + // Used in optimised AMG coarsening + virtual void expandAddrToZone(labelField&) const; + + //- Return the values of the given internal data adjacent to // the interface as a field virtual tmp interfaceInternalField diff --git a/src/foam/Make/files b/src/foam/Make/files index 4218eac70..f0f928e7d 100644 --- a/src/foam/Make/files +++ b/src/foam/Make/files @@ -141,13 +141,13 @@ $(StringStreams)/StringStreamsPrint.C Pstreams = $(Streams)/Pstreams $(Pstreams)/Pstream.C +$(Pstreams)/PstreamReduceOps.C $(Pstreams)/PstreamCommsStruct.C $(Pstreams)/PstreamGlobals.C $(Pstreams)/IPstream.C $(Pstreams)/OPstream.C $(Pstreams)/IPread.C $(Pstreams)/OPwrite.C -$(Pstreams)/PstreamsPrint.C dictionary = db/dictionary $(dictionary)/dictionary.C diff --git a/src/foam/containers/Lists/List/List.H b/src/foam/containers/Lists/List/List.H index ea96501d2..96119c4af 100644 --- a/src/foam/containers/Lists/List/List.H +++ b/src/foam/containers/Lists/List/List.H @@ -99,6 +99,7 @@ public: //- Return a null list inline static const List& null(); + // Constructors //- Null constructor. diff --git a/src/foam/db/IOstreams/Pstreams/IPread.C b/src/foam/db/IOstreams/Pstreams/IPread.C index d493491f2..abaddcef7 100644 --- a/src/foam/db/IOstreams/Pstreams/IPread.C +++ b/src/foam/db/IOstreams/Pstreams/IPread.C @@ -38,6 +38,8 @@ Foam::IPstream::IPstream const commsTypes commsType, const int fromProcNo, const label bufSize, + const int tag, + const label comm, streamFormat format, versionNumber version ) @@ -45,6 +47,8 @@ Foam::IPstream::IPstream Pstream(commsType, bufSize), Istream(format, version), fromProcNo_(fromProcNo), + tag_(tag), + comm_(comm), messageSize_(0) { setOpened(); @@ -52,17 +56,31 @@ Foam::IPstream::IPstream MPI_Status status; - // If the buffer size is not specified, probe the incomming message + // If the buffer size is not specified, probe the incoming message // and set it if (!bufSize) { - MPI_Probe(procID(fromProcNo_), msgType(), MPI_COMM_WORLD, &status); + MPI_Probe + ( + fromProcNo_, + tag_, + PstreamGlobals::MPICommunicators_[comm_], + &status + ); MPI_Get_count(&status, MPI_BYTE, &messageSize_); buf_.setSize(messageSize_); } - messageSize_ = read(commsType, fromProcNo_, buf_.begin(), buf_.size()); + messageSize_ = IPstream::read + ( + commsType, + fromProcNo_, + buf_.begin(), + buf_.size(), + tag_, + comm_ + ); if (!messageSize_) { @@ -83,9 +101,31 @@ Foam::label Foam::IPstream::read const commsTypes commsType, const int fromProcNo, char* buf, - const std::streamsize bufSize + const std::streamsize bufSize, + const int tag, + const label comm ) { + if (debug) + { + Pout<< "IPstream::read : starting read from:" << fromProcNo + << " tag:" << tag << " comm:" << comm + << " wanted size:" << label(bufSize) + << " commsType:" << Pstream::commsTypeNames[commsType] + << Foam::endl; + } + + if (Pstream::warnComm != -1 && comm != Pstream::warnComm) + { + Pout<< "IPstream::read : starting read from:" << fromProcNo + << " tag:" << tag << " comm:" << comm + << " wanted size:" << label(bufSize) + << " commsType:" << Pstream::commsTypeNames[commsType] + << " warnComm:" << Pstream::warnComm + << Foam::endl; + error::printStack(Pout); + } + if (commsType == blocking || commsType == scheduled) { MPI_Status status; @@ -96,10 +136,10 @@ Foam::label Foam::IPstream::read ( buf, bufSize, - MPI_PACKED, - procID(fromProcNo), - msgType(), - MPI_COMM_WORLD, + MPI_BYTE, + fromProcNo, + tag, + PstreamGlobals::MPICommunicators_[comm], &status ) ) @@ -144,10 +184,10 @@ Foam::label Foam::IPstream::read ( buf, bufSize, - MPI_PACKED, - procID(fromProcNo), - msgType(), - MPI_COMM_WORLD, + MPI_BYTE, + fromProcNo, + tag, + PstreamGlobals::MPICommunicators_[comm], &request ) ) @@ -162,8 +202,18 @@ Foam::label Foam::IPstream::read return 0; } - PstreamGlobals::IPstream_outstandingRequests_.append(request); + if (debug) + { + Pout<< "IPstream::read : started read from:" << fromProcNo + << " tag:" << tag << " read size:" << label(bufSize) + << " commsType:" << Pstream::commsTypeNames[commsType] + << " request:" << PstreamGlobals::outstandingRequests_.size() + << Foam::endl; + } + PstreamGlobals::outstandingRequests_.append(request); + + // Assume the message is completely received. return 1; } else @@ -180,56 +230,4 @@ Foam::label Foam::IPstream::read } -void Foam::IPstream::waitRequests() -{ - if (PstreamGlobals::IPstream_outstandingRequests_.size()) - { - if - ( - MPI_Waitall - ( - PstreamGlobals::IPstream_outstandingRequests_.size(), - PstreamGlobals::IPstream_outstandingRequests_.begin(), - MPI_STATUSES_IGNORE - ) - ) - { - FatalErrorIn - ( - "IPstream::waitRequests()" - ) << "MPI_Waitall returned with error" << endl; - } - - PstreamGlobals::IPstream_outstandingRequests_.clear(); - } -} - - -bool Foam::IPstream::finishedRequest(const label i) -{ - if (i >= PstreamGlobals::IPstream_outstandingRequests_.size()) - { - FatalErrorIn - ( - "IPstream::finishedRequest(const label)" - ) << "There are " - << PstreamGlobals::IPstream_outstandingRequests_.size() - << " outstanding send requests and you are asking for i=" << i - << nl - << "Maybe you are mixing blocking/non-blocking comms?" - << Foam::abort(FatalError); - } - - int flag; - MPI_Test - ( - &PstreamGlobals::IPstream_outstandingRequests_[i], - &flag, - MPI_STATUS_IGNORE - ); - - return flag != 0; -} - - // ************************************************************************* // diff --git a/src/foam/db/IOstreams/Pstreams/IPreadToken.C b/src/foam/db/IOstreams/Pstreams/IPreadToken.C deleted file mode 100644 index d4b8d228c..000000000 --- a/src/foam/db/IOstreams/Pstreams/IPreadToken.C +++ /dev/null @@ -1,187 +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 . - -\*---------------------------------------------------------------------------*/ - -#include "IPstream.H" -#include "token.H" -#include - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -Istream& IPstream::read(token& t) -{ - // Return the put back token if it exists - if (Istream::getBack(t)) - { - return *this; - } - - char c; - - // return on error - if (!read(c)) - { - t.setBad(); - return *this; - } - - // Set the line number of this token to the current stream line number - t.lineNumber() = lineNumber(); - - // Analyse input starting with this character. - switch (c) - { - // Punctuation - case token::END_STATEMENT : - case token::BEGIN_LIST : - case token::END_LIST : - case token::BEGIN_SQR : - case token::END_SQR : - case token::BEGIN_BLOCK : - case token::END_BLOCK : - case token::COLON : - case token::COMMA : - case token::ASSIGN : - case token::ADD : - case token::SUBTRACT : - case token::MULTIPLY : - case token::DIVIDE : - { - t = token::punctuationToken(c); - return *this; - } - - // Word - case token::WORD : - { - word* wPtr = new word; - if (read(*wPtr)) - { - if (token::compound::isCompound(*wPtr)) - { - t = token::compound::New(*wPtr, *this).ptr(); - delete wPtr; - } - else - { - t = wPtr; - } - } - else - { - delete wPtr; - t.setBad(); - } - return *this; - } - - // String - case token::STRING : - { - string* sPtr = new string; - if (read(*sPtr)) - { - t = sPtr; - } - else - { - delete sPtr; - t.setBad(); - } - return *this; - } - - // Label - case token::LABEL : - { - label l; - if (read(l)) - { - t = l; - } - else - { - t.setBad(); - } - return *this; - } - - // floatScalar - case token::FLOAT_SCALAR : - { - floatScalar s; - if (read(s)) - { - t = s; - } - else - { - t.setBad(); - } - return *this; - } - - // doubleScalar - case token::DOUBLE_SCALAR : - { - doubleScalar s; - if (read(s)) - { - t = s; - } - else - { - t.setBad(); - } - return *this; - } - - // Character (returned as a single character word) or error - default: - { - if (isalpha(c)) - { - t = word(c); - return *this; - } - - setBad(); - t.setBad(); - - return *this; - } - } -} - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// ************************************************************************* // diff --git a/src/foam/db/IOstreams/Pstreams/IPstream.C b/src/foam/db/IOstreams/Pstreams/IPstream.C index b3b31946d..090c1ceb6 100644 --- a/src/foam/db/IOstreams/Pstreams/IPstream.C +++ b/src/foam/db/IOstreams/Pstreams/IPstream.C @@ -324,4 +324,13 @@ Foam::Istream& Foam::IPstream::rewind() } +void Foam::IPstream::print(Ostream& os) const +{ + os << "Reading from processor " << fromProcNo_ + << " using communicator " << comm_ + << " and tag " << tag_ + << Foam::endl; +} + + // ************************************************************************* // diff --git a/src/foam/db/IOstreams/Pstreams/IPstream.H b/src/foam/db/IOstreams/Pstreams/IPstream.H index c69b024bc..bac005a20 100644 --- a/src/foam/db/IOstreams/Pstreams/IPstream.H +++ b/src/foam/db/IOstreams/Pstreams/IPstream.H @@ -58,6 +58,12 @@ class IPstream //- ID of sending processor int fromProcNo_; + //- Message tag + const int tag_; + + //- Communicator + const label comm_; + //- Message size label messageSize_; @@ -86,12 +92,14 @@ public: const commsTypes commsType, const int fromProcNo, const label bufSize = 0, + const int tag = Pstream::msgType(), + const label comm = Pstream::worldComm, streamFormat format = BINARY, versionNumber version = currentVersion ); - // Destructor + //- Destructor ~IPstream(); @@ -115,15 +123,11 @@ public: const commsTypes commsType, const int fromProcNo, char* buf, - const std::streamsize bufSize + const std::streamsize bufSize, + const int tag = Pstream::msgType(), + const label communicator = 0 ); - //- Non-blocking receives: wait until all have finished. - static void waitRequests(); - - //- Non-blocking receives: has request i finished? - static bool finishedRequest(const label i); - //- Return next token from stream Istream& read(token&); diff --git a/src/foam/db/IOstreams/Pstreams/OPstream.C b/src/foam/db/IOstreams/Pstreams/OPstream.C index 204330e6b..d643f27f1 100644 --- a/src/foam/db/IOstreams/Pstreams/OPstream.C +++ b/src/foam/db/IOstreams/Pstreams/OPstream.C @@ -92,13 +92,17 @@ Foam::OPstream::OPstream const commsTypes commsType, const int toProcNo, const label bufSize, + const int tag, + const label comm, streamFormat format, versionNumber version ) : Pstream(commsType, bufSize), Ostream(format, version), - toProcNo_(toProcNo) + toProcNo_(toProcNo), + tag_(tag), + comm_(comm) { setOpened(); setGood(); @@ -233,4 +237,12 @@ Foam::Ostream& Foam::OPstream::write(const char* data, std::streamsize count) } +void Foam::OPstream::print(Ostream& os) const +{ + os << "Writing from processor " << toProcNo_ + << " to processor " << myProcNo() << " in communicator " << comm_ + << " and tag " << tag_ << Foam::endl; +} + + // ************************************************************************* // diff --git a/src/foam/db/IOstreams/Pstreams/OPstream.H b/src/foam/db/IOstreams/Pstreams/OPstream.H index 623637618..c68e563f8 100644 --- a/src/foam/db/IOstreams/Pstreams/OPstream.H +++ b/src/foam/db/IOstreams/Pstreams/OPstream.H @@ -58,6 +58,12 @@ class OPstream // ID of receiving processor int toProcNo_; + //- Message tag + const int tag_; + + //- Communicator + const label comm_; + // Private member functions @@ -88,14 +94,15 @@ public: const commsTypes commsType, const int toProcNo, const label bufSize = 0, + const int tag = Pstream::msgType(), + const label comm = Pstream::worldComm, streamFormat format = BINARY, versionNumber version = currentVersion ); - // Destructor - - ~OPstream(); + //- Destructor + ~OPstream(); // Member functions @@ -117,15 +124,11 @@ public: const commsTypes commsType, const int toProcNo, const char* buf, - const std::streamsize bufSize + const std::streamsize bufSize, + const int tag = Pstream::msgType(), + const label comm = Pstream::worldComm ); - //- Non-blocking writes: wait until all have finished. - static void waitRequests(); - - //- Non-blocking writes: has request i finished? - static bool finishedRequest(const label i); - //- Write next token to stream Ostream& write(const token&); diff --git a/src/foam/db/IOstreams/Pstreams/OPwrite.C b/src/foam/db/IOstreams/Pstreams/OPwrite.C index 644a83b00..5f3302bc9 100644 --- a/src/foam/db/IOstreams/Pstreams/OPwrite.C +++ b/src/foam/db/IOstreams/Pstreams/OPwrite.C @@ -42,7 +42,9 @@ Foam::OPstream::~OPstream() commsType_, toProcNo_, buf_.begin(), - bufPosition_ + bufPosition_, + tag_, + comm_ ) ) { @@ -60,9 +62,34 @@ bool Foam::OPstream::write const commsTypes commsType, const int toProcNo, const char* buf, - const std::streamsize bufSize + const std::streamsize bufSize, + const int tag, + const label comm ) { + if (debug) + { + Pout<< "OPstream::write : starting write to:" << toProcNo + << " tag:" << tag + << " comm:" << comm << " size:" << label(bufSize) + << " commsType:" << Pstream::commsTypeNames[commsType] + << Foam::endl; + } + if (Pstream::warnComm != -1 && comm != Pstream::warnComm) + { + Pout<< "OPstream::write : starting write to:" << toProcNo + << " tag:" << tag + << " comm:" << comm << " size:" << label(bufSize) + << " commsType:" << Pstream::commsTypeNames[commsType] + << " warnComm:" << Pstream::warnComm + << Foam::endl; + + error::printStack(Pout); + } + + + PstreamGlobals::checkCommunicator(comm, toProcNo); + bool transferFailed = true; if (commsType == blocking) @@ -71,11 +98,19 @@ bool Foam::OPstream::write ( const_cast(buf), bufSize, - MPI_PACKED, - procID(toProcNo), - msgType(), - MPI_COMM_WORLD + MPI_BYTE, + toProcNo, //procID(toProcNo), + tag, + PstreamGlobals::MPICommunicators_[comm] // MPI_COMM_WORLD ); + + if (debug) + { + Pout<< "OPstream::write : finished write to:" << toProcNo + << " tag:" << tag << " size:" << label(bufSize) + << " commsType:" << Pstream::commsTypeNames[commsType] + << Foam::endl; + } } else if (commsType == scheduled) { @@ -83,11 +118,19 @@ bool Foam::OPstream::write ( const_cast(buf), bufSize, - MPI_PACKED, - procID(toProcNo), - msgType(), - MPI_COMM_WORLD + MPI_BYTE, + toProcNo, //procID(toProcNo), + tag, + PstreamGlobals::MPICommunicators_[comm] // MPI_COMM_WORLD ); + + if (debug) + { + Pout<< "OPstream::write : finished write to:" << toProcNo + << " tag:" << tag << " size:" << label(bufSize) + << " commsType:" << Pstream::commsTypeNames[commsType] + << Foam::endl; + } } else if (commsType == nonBlocking) { @@ -97,14 +140,23 @@ bool Foam::OPstream::write ( const_cast(buf), bufSize, - MPI_PACKED, - procID(toProcNo), - msgType(), - MPI_COMM_WORLD, + MPI_BYTE, + toProcNo, //procID(toProcNo), + tag, + PstreamGlobals::MPICommunicators_[comm],// MPI_COMM_WORLD, &request ); - PstreamGlobals::OPstream_outstandingRequests_.append(request); + if (debug) + { + Pout<< "OPstream::write : started write to:" << toProcNo + << " tag:" << tag << " size:" << label(bufSize) + << " commsType:" << Pstream::commsTypeNames[commsType] + << " request:" << PstreamGlobals::outstandingRequests_.size() + << Foam::endl; + } + + PstreamGlobals::outstandingRequests_.append(request); } else { @@ -120,56 +172,4 @@ bool Foam::OPstream::write } -void Foam::OPstream::waitRequests() -{ - if (PstreamGlobals::OPstream_outstandingRequests_.size()) - { - if - ( - MPI_Waitall - ( - PstreamGlobals::OPstream_outstandingRequests_.size(), - PstreamGlobals::OPstream_outstandingRequests_.begin(), - MPI_STATUSES_IGNORE - ) - ) - { - FatalErrorIn - ( - "OPstream::waitRequests()" - ) << "MPI_Waitall returned with error" << Foam::endl; - } - - PstreamGlobals::OPstream_outstandingRequests_.clear(); - } -} - - -bool Foam::OPstream::finishedRequest(const label i) -{ - if (i >= PstreamGlobals::OPstream_outstandingRequests_.size()) - { - FatalErrorIn - ( - "OPstream::finishedRequest(const label)" - ) << "There are " - << PstreamGlobals::OPstream_outstandingRequests_.size() - << " outstanding send requests and you are asking for i=" << i - << nl - << "Maybe you are mixing blocking/non-blocking comms?" - << Foam::abort(FatalError); - } - - int flag; - MPI_Test - ( - &PstreamGlobals::OPstream_outstandingRequests_[i], - &flag, - MPI_STATUS_IGNORE - ); - - return flag != 0; -} - - // ************************************************************************* // diff --git a/src/foam/db/IOstreams/Pstreams/Pstream.C b/src/foam/db/IOstreams/Pstreams/Pstream.C index 422aff9c2..a7e6dc63c 100644 --- a/src/foam/db/IOstreams/Pstreams/Pstream.C +++ b/src/foam/db/IOstreams/Pstreams/Pstream.C @@ -23,11 +23,6 @@ License \*---------------------------------------------------------------------------*/ - -#include -#include -#include - #include "mpi.h" #include "Pstream.H" @@ -35,15 +30,12 @@ License #include "debug.H" #include "dictionary.H" #include "OSspecific.H" +#include "PstreamGlobals.H" +#include "SubList.H" -#if defined(WM_SP) -# define MPI_SCALAR MPI_FLOAT -#elif defined(WM_DP) -# define MPI_SCALAR MPI_DOUBLE -#elif defined(WM_LDP) -# define MPI_SCALAR MPI_LONG_DOUBLE -#endif - +#include +#include +#include // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -63,65 +55,85 @@ const Foam::NamedEnum // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::Pstream::setParRun() +void Foam::Pstream::setParRun(const label nProcs) { - parRun_ = true; + if (nProcs == 0) + { + parRun_ = false; + freeCommunicator(Pstream::worldComm); - Pout.prefix() = '[' + name(myProcNo()) + "] "; - Perr.prefix() = '[' + name(myProcNo()) + "] "; + label comm = allocateCommunicator(-1, labelList(1, label(0)), false); + + if (comm != Pstream::worldComm) + { + FatalErrorIn("Pstream::setParRun(const label)") + << "problem : comm:" << comm + << " Pstream::worldComm:" << Pstream::worldComm + << Foam::exit(FatalError); + } + + Pout.prefix() = ""; + Perr.prefix() = ""; + } + else + { + parRun_ = true; + + // Redo worldComm communicator (created at static initialisation) + freeCommunicator(Pstream::worldComm); + label comm = allocateCommunicator(-1, identity(nProcs), true); + + if (comm != Pstream::worldComm) + { + FatalErrorIn("Pstream::setParRun(const label)") + << "problem : comm:" << comm + << " Pstream::worldComm:" << Pstream::worldComm + << Foam::exit(FatalError); + } + + Pout.prefix() = '[' + name(myProcNo()) + "] "; + Perr.prefix() = '[' + name(myProcNo()) + "] "; + } } -void Foam::Pstream::calcLinearComm(const label nProcs) +Foam::List Foam::Pstream::calcLinearComm +( + const label nProcs +) { - linearCommunication_.setSize(nProcs); + List linearCommunication(nProcs); // Master labelList belowIDs(nProcs - 1); - forAll(belowIDs, i) + forAll (belowIDs, i) { belowIDs[i] = i + 1; } - linearCommunication_[0] = commsStruct + linearCommunication[0] = commsStruct ( nProcs, 0, -1, belowIDs, - labelList(0) + labelList() ); // Slaves. Have no below processors, only communicate up to master for (label procID = 1; procID < nProcs; procID++) { - linearCommunication_[procID] = commsStruct + linearCommunication[procID] = commsStruct ( nProcs, procID, 0, - labelList(0), - labelList(0) + labelList(), + labelList() ); } -} - -// Append my children (and my children children etc.) to allReceives. -void Foam::Pstream::collectReceives -( - const label procID, - const List& receives, - dynamicLabelList& allReceives -) -{ - const dynamicLabelList& myChildren = receives[procID]; - - forAll(myChildren, childI) - { - allReceives.append(myChildren[childI]); - collectReceives(myChildren[childI], receives, allReceives); - } + return linearCommunication; } @@ -151,7 +163,10 @@ void Foam::Pstream::collectReceives // 5 - 4 // 6 7 4 // 7 - 6 -void Foam::Pstream::calcTreeComm(label nProcs) +Foam::List Foam::Pstream::calcTreeComm +( + const label nProcs +) { label nLevels = 1; while ((1 << nLevels) < nProcs) @@ -159,11 +174,9 @@ void Foam::Pstream::calcTreeComm(label nProcs) nLevels++; } - List receives(nProcs); + List receives(nProcs); labelList sends(nProcs, -1); - // Info<< "Using " << nLevels << " communication levels" << endl; - label offset = 2; label childOffset = offset/2; @@ -190,18 +203,18 @@ void Foam::Pstream::calcTreeComm(label nProcs) // For all processors find the processors it receives data from // (and the processors they receive data from etc.) - List allReceives(nProcs); + List allReceives(nProcs); for (label procID = 0; procID < nProcs; procID++) { collectReceives(procID, receives, allReceives[procID]); } - treeCommunication_.setSize(nProcs); + List treeCommunication(nProcs); for (label procID = 0; procID < nProcs; procID++) { - treeCommunication_[procID] = commsStruct + treeCommunication[procID] = commsStruct ( nProcs, procID, @@ -210,6 +223,26 @@ void Foam::Pstream::calcTreeComm(label nProcs) allReceives[procID].shrink() ); } + + return treeCommunication; +} + + +// Append my children (and my children;s children etc.) to allReceives. +void Foam::Pstream::collectReceives +( + const label procID, + const List& receives, + dynamicLabelList& allReceives +) +{ + const dynamicLabelList& myChildren = receives[procID]; + + forAll (myChildren, childI) + { + allReceives.append(myChildren[childI]); + collectReceives(myChildren[childI], receives, allReceives); + } } @@ -222,6 +255,292 @@ void Foam::Pstream::initCommunicationSchedule() } +void Foam::Pstream::allocatePstreamCommunicator +( + const label parentIndex, + const label index +) +{ + if (index == PstreamGlobals::MPIGroups_.size()) + { + // Extend storage with dummy values + MPI_Group newGroup; + PstreamGlobals::MPIGroups_.append(newGroup); + MPI_Comm newComm; + PstreamGlobals::MPICommunicators_.append(newComm); + } + else if (index > PstreamGlobals::MPIGroups_.size()) + { + FatalErrorIn + ( + "Pstream::allocatePstreamCommunicator\n" + "(\n" + " const label parentIndex,\n" + " const labelList& subRanks\n" + ")\n" + ) << "PstreamGlobals out of sync with Pstream data. Problem." + << Foam::exit(FatalError); + } + + + if (parentIndex == -1) + { + // Allocate world communicator + + if (index != Pstream::worldComm) + { + FatalErrorIn + ( + "Pstream::allocatePstreamCommunicator\n" + "(\n" + " const label parentIndex,\n" + " const labelList& subRanks\n" + ")\n" + ) << "world communicator should always be index " + << Pstream::worldComm << Foam::exit(FatalError); + } + + PstreamGlobals::MPICommunicators_[index] = MPI_COMM_WORLD; + MPI_Comm_group(MPI_COMM_WORLD, &PstreamGlobals::MPIGroups_[index]); + MPI_Comm_rank + ( + PstreamGlobals::MPICommunicators_[index], + &myProcNo_[index] + ); + + // Set the number of processes to the actual number + int numProcs; + MPI_Comm_size(PstreamGlobals::MPICommunicators_[index], &numProcs); + procIDs_[index] = identity(numProcs); + } + else + { + // Create new group + MPI_Group_incl + ( + PstreamGlobals::MPIGroups_[parentIndex], + procIDs_[index].size(), + procIDs_[index].begin(), + &PstreamGlobals::MPIGroups_[index] + ); + + // Create new communicator + MPI_Comm_create + ( + PstreamGlobals::MPICommunicators_[parentIndex], + PstreamGlobals::MPIGroups_[index], + &PstreamGlobals::MPICommunicators_[index] + ); + + if (PstreamGlobals::MPICommunicators_[index] == MPI_COMM_NULL) + { + myProcNo_[index] = -1; + } + else + { + if + ( + MPI_Comm_rank + ( + PstreamGlobals::MPICommunicators_[index], + &myProcNo_[index] + ) + ) + { + FatalErrorIn + ( + "Pstream::allocatePstreamCommunicator\n" + "(\n" + " const label,\n" + " const labelList&\n" + ")\n" + ) << "Problem :" + << " when allocating communicator at " << index + << " from ranks " << procIDs_[index] + << " of parent " << parentIndex + << " cannot find my own rank" + << Foam::exit(FatalError); + } + } + } +} + + +void Foam::Pstream::freePstreamCommunicator(const label communicator) +{ + if (communicator != Pstream::worldComm) + { + if (PstreamGlobals::MPICommunicators_[communicator] != MPI_COMM_NULL) + { + // Free communicator. Sets communicator to MPI_COMM_NULL + MPI_Comm_free(&PstreamGlobals::MPICommunicators_[communicator]); + } + if (PstreamGlobals::MPIGroups_[communicator] != MPI_GROUP_NULL) + { + // Free greoup. Sets group to MPI_GROUP_NULL + MPI_Group_free(&PstreamGlobals::MPIGroups_[communicator]); + } + } +} + + +Foam::label Foam::Pstream::allocateCommunicator +( + const label parentIndex, + const labelList& subRanks, + const bool doPstream +) +{ + label index; + if (!freeComms_.empty()) + { + index = freeComms_.pop(); + } + else + { + // Extend storage + index = parentCommunicator_.size(); + + myProcNo_.append(-1); + procIDs_.append(List()); + parentCommunicator_.append(-1); + linearCommunication_.append(List()); + treeCommunication_.append(List()); + } + + if (debug) + { + Pout<< "Communicators : Allocating communicator " << index << endl + << " parent : " << parentIndex << endl + << " procs : " << subRanks << endl + << endl; + } + + // Initialise; overwritten by allocatePstreamCommunicator + myProcNo_[index] = 0; + + // Convert from label to int + procIDs_[index].setSize(subRanks.size()); + forAll (procIDs_[index], i) + { + procIDs_[index][i] = subRanks[i]; + + // Enforce incremental order (so index is rank in next communicator) + if (i >= 1 && subRanks[i] <= subRanks[i - 1]) + { + FatalErrorIn + ( + "Pstream::allocateCommunicator" + "(const label, const labelList&, const bool)" + ) << "subranks not sorted : " << subRanks + << " when allocating subcommunicator from parent " + << parentIndex + << Foam::abort(FatalError); + } + } + parentCommunicator_[index] = parentIndex; + + linearCommunication_[index] = calcLinearComm(procIDs_[index].size()); + treeCommunication_[index] = calcTreeComm(procIDs_[index].size()); + + + if (doPstream && parRun()) + { + allocatePstreamCommunicator(parentIndex, index); + } + + return index; +} + + +void Foam::Pstream::freeCommunicator +( + const label communicator, + const bool doPstream +) +{ + if (debug) + { + Pout<< "Communicators : Freeing communicator " << communicator << endl + << " parent : " << parentCommunicator_[communicator] << endl + << " myProcNo : " << myProcNo_[communicator] << endl + << endl; + } + + if (doPstream && parRun()) + { + freePstreamCommunicator(communicator); + } + myProcNo_[communicator] = -1; + //procIDs_[communicator].clear(); + parentCommunicator_[communicator] = -1; + linearCommunication_[communicator].clear(); + treeCommunication_[communicator].clear(); + + freeComms_.push(communicator); +} + + +void Foam::Pstream::freeCommunicators(const bool doPstream) +{ + forAll (myProcNo_, communicator) + { + if (myProcNo_[communicator] != -1) + { + freeCommunicator(communicator, doPstream); + } + } +} + + +int Foam::Pstream::baseProcNo(const label myComm, const int myProcID) +{ + int procID = myProcID; + label comm = myComm; + + while (parent(comm) != -1) + { + const List& parentRanks = Pstream::procID(comm); + procID = parentRanks[procID]; + comm = Pstream::parent(comm); + } + + return procID; +} + + +Foam::label Foam::Pstream::procNo(const label myComm, const int baseProcID) +{ + const List& parentRanks = procID(myComm); + label parentComm = parent(myComm); + + if (parentComm == -1) + { + return findIndex(parentRanks, baseProcID); + } + else + { + label parentRank = procNo(parentComm, baseProcID); + return findIndex(parentRanks, parentRank); + } +} + + +Foam::label Foam::Pstream::procNo +( + const label myComm, + const label currentComm, + const int currentProcID +) +{ + label physProcID = Pstream::baseProcNo(currentComm, currentProcID); + return procNo(myComm, physProcID); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + + // NOTE: // valid parallel options vary between implementations, but flag common ones. // if they are not removed by MPI_Init(), the subsequent argument processing @@ -233,7 +552,6 @@ void Foam::Pstream::addValidParOptions(HashTable& validParOptions) validParOptions.insert("p4wd", "directory"); validParOptions.insert("p4amslave", ""); validParOptions.insert("p4yourname", "hostname"); - validParOptions.insert("GAMMANP", "number of instances"); validParOptions.insert("machinefile", "machine file"); } @@ -244,24 +562,25 @@ bool Foam::Pstream::init(int& argc, char**& argv) int numprocs; MPI_Comm_size(MPI_COMM_WORLD, &numprocs); - MPI_Comm_rank(MPI_COMM_WORLD, &myProcNo_); + int myRank; + MPI_Comm_rank(MPI_COMM_WORLD, &myRank); + + if (debug) + { + Pout<< "Pstream::init : initialised with numProcs:" << numprocs + << " myRank:" << myRank << endl; + } if (numprocs <= 1) { FatalErrorIn("Pstream::init(int& argc, char**& argv)") - << "bool Pstream::init(int& argc, char**& argv) : " + << "bool IPstream::init(int& argc, char**& argv) : " "attempt to run parallel on 1 processor" << Foam::abort(FatalError); } - procIDs_.setSize(numprocs); - - forAll(procIDs_, procNo) - { - procIDs_[procNo] = procNo; - } - - setParRun(); + // Initialise parallel structure + setParRun(numprocs); # ifndef SGIMPI string bufferSizeName = getEnv("MPI_BUFFER_SIZE"); @@ -284,15 +603,12 @@ bool Foam::Pstream::init(int& argc, char**& argv) } # endif - int processorNameLen; - char processorName[MPI_MAX_PROCESSOR_NAME]; - - MPI_Get_processor_name(processorName, &processorNameLen); - - //signal(SIGABRT, stop); - - // Now that nprocs is known construct communication tables. - initCommunicationSchedule(); + //int processorNameLen; + //char processorName[MPI_MAX_PROCESSOR_NAME]; + // + //MPI_Get_processor_name(processorName, &processorNameLen); + //processorName[processorNameLen] = '\0'; + //Pout<< "Processor name:" << processorName << endl; return true; } @@ -300,6 +616,11 @@ bool Foam::Pstream::init(int& argc, char**& argv) void Foam::Pstream::exit(int errnum) { + if (debug) + { + Pout<< "Pstream::exit." << endl; + } + # ifndef SGIMPI int size; char* buff; @@ -307,6 +628,28 @@ void Foam::Pstream::exit(int errnum) delete[] buff; # endif + if (PstreamGlobals::outstandingRequests_.size()) + { + label n = PstreamGlobals::outstandingRequests_.size(); + PstreamGlobals::outstandingRequests_.clear(); + + WarningIn("Pstream::exit(int)") + << "There are still " << n << " outstanding MPI_Requests." << endl + << "This means that your code exited before doing a" + << " Pstream::waitRequests()." << endl + << "This should not happen for a normal code exit." + << endl; + } + + // Clean mpi communicators + forAll (myProcNo_, communicator) + { + if (myProcNo_[communicator] != -1) + { + freePstreamCommunicator(communicator); + } + } + if (errnum == 0) { MPI_Finalize(); @@ -325,316 +668,244 @@ void Foam::Pstream::abort() } -void Foam::reduce(scalar& Value, const sumOp& bop) +Foam::label Foam::Pstream::nRequests() { - if (!Pstream::parRun()) + return PstreamGlobals::outstandingRequests_.size(); +} + + +void Foam::Pstream::resetRequests(const label i) +{ + if (i < PstreamGlobals::outstandingRequests_.size()) { - return; + PstreamGlobals::outstandingRequests_.setSize(i); + } +} + + +void Foam::Pstream::waitRequests(const label start) +{ + if (debug) + { + Pout<< "Pstream::waitRequests : starting wait for " + << PstreamGlobals::outstandingRequests_.size()-start + << " outstanding requests starting at " << start << endl; } - if (Pstream::nProcs() <= Pstream::nProcsSimpleSum()) + if (PstreamGlobals::outstandingRequests_.size()) { - if (Pstream::master()) - { - for - ( - int slave=Pstream::firstSlave(); - slave<=Pstream::lastSlave(); - slave++ - ) - { - scalar value; + SubList waitRequests + ( + PstreamGlobals::outstandingRequests_, + PstreamGlobals::outstandingRequests_.size() - start, + start + ); - if - ( - MPI_Recv - ( - &value, - 1, - MPI_SCALAR, - Pstream::procID(slave), - Pstream::msgType(), - MPI_COMM_WORLD, - MPI_STATUS_IGNORE - ) - ) - { - FatalErrorIn - ( - "reduce(scalar& Value, const sumOp& sumOp)" - ) << "MPI_Recv failed" - << Foam::abort(FatalError); - } - - Value = bop(Value, value); - } - } - else - { - if + if + ( + MPI_Waitall ( - MPI_Send - ( - &Value, - 1, - MPI_SCALAR, - Pstream::procID(Pstream::masterNo()), - Pstream::msgType(), - MPI_COMM_WORLD - ) + waitRequests.size(), + waitRequests.begin(), + MPI_STATUSES_IGNORE ) - { - FatalErrorIn - ( - "reduce(scalar& Value, const sumOp& sumOp)" - ) << "MPI_Send failed" - << Foam::abort(FatalError); - } + ) + { + FatalErrorIn + ( + "Pstream::waitRequests()" + ) << "MPI_Waitall returned with error" << Foam::endl; } + resetRequests(start); + } - if (Pstream::master()) - { - for - ( - int slave=Pstream::firstSlave(); - slave<=Pstream::lastSlave(); - slave++ - ) - { - if - ( - MPI_Send - ( - &Value, - 1, - MPI_SCALAR, - Pstream::procID(slave), - Pstream::msgType(), - MPI_COMM_WORLD - ) - ) - { - FatalErrorIn - ( - "reduce(scalar& Value, const sumOp& sumOp)" - ) << "MPI_Send failed" - << Foam::abort(FatalError); - } - } - } - else - { - if - ( - MPI_Recv - ( - &Value, - 1, - MPI_SCALAR, - Pstream::procID(Pstream::masterNo()), - Pstream::msgType(), - MPI_COMM_WORLD, - MPI_STATUS_IGNORE - ) - ) - { - FatalErrorIn - ( - "reduce(scalar& Value, const sumOp& sumOp)" - ) << "MPI_Recv failed" - << Foam::abort(FatalError); - } - } + if (debug) + { + Pout<< "Pstream::waitRequests : finished wait." << endl; + } +} + + +void Foam::Pstream::waitRequest(const label i) +{ + if (debug) + { + Pout<< "Pstream::waitRequest : starting wait for request:" << i + << endl; + } + + if (i >= PstreamGlobals::outstandingRequests_.size()) + { + FatalErrorIn + ( + "Pstream::waitRequest(const label)" + ) << "There are " << PstreamGlobals::outstandingRequests_.size() + << " outstanding send requests and you are asking for i=" << i + << nl + << "Maybe you are mixing blocking/non-blocking comms?" + << Foam::abort(FatalError); + } + + if + ( + MPI_Wait + ( + &PstreamGlobals::outstandingRequests_[i], + MPI_STATUS_IGNORE + ) + ) + { + FatalErrorIn + ( + "Pstream::waitRequest()" + ) << "MPI_Wait returned with error" << Foam::endl; + } + + if (debug) + { + Pout<< "Pstream::waitRequest : finished wait for request:" << i + << endl; + } +} + + +bool Foam::Pstream::finishedRequest(const label i) +{ + if (debug) + { + Pout<< "Pstream::finishedRequest : checking request:" << i + << endl; + } + + if (i >= PstreamGlobals::outstandingRequests_.size()) + { + FatalErrorIn + ( + "Pstream::finishedRequest(const label)" + ) << "There are " << PstreamGlobals::outstandingRequests_.size() + << " outstanding send requests and you are asking for i=" << i + << nl + << "Maybe you are mixing blocking/non-blocking comms?" + << Foam::abort(FatalError); + } + + int flag; + MPI_Test + ( + &PstreamGlobals::outstandingRequests_[i], + &flag, + MPI_STATUS_IGNORE + ); + + if (debug) + { + Pout<< "Pstream::finishedRequest : finished request:" << i + << endl; + } + + return flag != 0; +} + + +int Foam::Pstream::allocateTag(const char* s) +{ + int tag; + if (PstreamGlobals::freedTags_.size()) + { + tag = PstreamGlobals::freedTags_.remove(); } else { - scalar sum; - MPI_Allreduce(&Value, &sum, 1, MPI_SCALAR, MPI_SUM, MPI_COMM_WORLD); - Value = sum; - - /* - int myProcNo = Pstream::myProcNo(); - int nProcs = Pstream::nProcs(); - - // - // receive from children - // - int level = 1; - int thisLevelOffset = 2; - int childLevelOffset = thisLevelOffset/2; - int childProcId = 0; - - while - ( - (childLevelOffset < nProcs) - && (myProcNo % thisLevelOffset) == 0 - ) - { - childProcId = myProcNo + childLevelOffset; - - scalar value; - - if (childProcId < nProcs) - { - if - ( - MPI_Recv - ( - &value, - 1, - MPI_SCALAR, - Pstream::procID(childProcId), - Pstream::msgType(), - MPI_COMM_WORLD, - MPI_STATUS_IGNORE - ) - ) - { - FatalErrorIn - ( - "reduce(scalar& Value, const sumOp& sumOp)" - ) << "MPI_Recv failed" - << Foam::abort(FatalError); - } - - Value = bop(Value, value); - } - - level++; - thisLevelOffset <<= 1; - childLevelOffset = thisLevelOffset/2; - } - - // - // send and receive from parent - // - if (!Pstream::master()) - { - int parentId = myProcNo - (myProcNo % thisLevelOffset); - - if - ( - MPI_Send - ( - &Value, - 1, - MPI_SCALAR, - Pstream::procID(parentId), - Pstream::msgType(), - MPI_COMM_WORLD - ) - ) - { - FatalErrorIn - ( - "reduce(scalar& Value, const sumOp& sumOp)" - ) << "MPI_Send failed" - << Foam::abort(FatalError); - } - - if - ( - MPI_Recv - ( - &Value, - 1, - MPI_SCALAR, - Pstream::procID(parentId), - Pstream::msgType(), - MPI_COMM_WORLD, - MPI_STATUS_IGNORE - ) - ) - { - FatalErrorIn - ( - "reduce(scalar& Value, const sumOp& sumOp)" - ) << "MPI_Recv failed" - << Foam::abort(FatalError); - } - } - - - // - // distribute to my children - // - level--; - thisLevelOffset >>= 1; - childLevelOffset = thisLevelOffset/2; - - while (level > 0) - { - childProcId = myProcNo + childLevelOffset; - - if (childProcId < nProcs) - { - if - ( - MPI_Send - ( - &Value, - 1, - MPI_SCALAR, - Pstream::procID(childProcId), - Pstream::msgType(), - MPI_COMM_WORLD - ) - ) - { - FatalErrorIn - ( - "reduce(scalar& Value, const sumOp& sumOp)" - ) << "MPI_Send failed" - << Foam::abort(FatalError); - } - } - - level--; - thisLevelOffset >>= 1; - childLevelOffset = thisLevelOffset/2; - } - */ + tag = PstreamGlobals::nTags_++; } + + if (debug) + { + Pout<< "Pstream::allocateTag " << s + << " : tag:" << tag + << endl; + } + + return tag; +} + + +int Foam::Pstream::allocateTag(const word& s) +{ + int tag; + if (PstreamGlobals::freedTags_.size()) + { + tag = PstreamGlobals::freedTags_.remove(); + } + else + { + tag = PstreamGlobals::nTags_++; + } + + if (debug) + { + Pout<< "Pstream::allocateTag " << s + << " : tag:" << tag + << endl; + } + + return tag; +} + + +void Foam::Pstream::freeTag(const char* s, const int tag) +{ + if (debug) + { + Pout<< "Pstream::freeTag " << s << " tag:" << tag << endl; + } + PstreamGlobals::freedTags_.append(tag); +} + + +void Foam::Pstream::freeTag(const word& s, const int tag) +{ + if (debug) + { + Pout<< "Pstream::freeTag " << s << " tag:" << tag << endl; + } + PstreamGlobals::freedTags_.append(tag); } // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -// Initialise my process number to 0 (the master) -int Foam::Pstream::myProcNo_(0); - - // By default this is not a parallel run bool Foam::Pstream::parRun_(false); +// Free communicators +Foam::LIFOStack Foam::Pstream::freeComms_; + +// My processor number +Foam::DynamicList Foam::Pstream::myProcNo_(10); // List of process IDs -Foam::List Foam::Pstream::procIDs_(1, 0); +Foam::DynamicList > Foam::Pstream::procIDs_(10); +// Parent communicator +Foam::DynamicList Foam::Pstream::parentCommunicator_(10); // Standard transfer message type const int Foam::Pstream::msgType_(1); - // Linear communication schedule -Foam::List Foam::Pstream::linearCommunication_(0); - +Foam::DynamicList > +Foam::Pstream::linearCommunication_(10); // Multi level communication schedule -Foam::List Foam::Pstream::treeCommunication_(0); +Foam::DynamicList > +Foam::Pstream::treeCommunication_(10); -// Should compact transfer be used in which floats replace doubles -// reducing the bandwidth requirement at the expense of some loss -// in accuracy -const Foam::debug::optimisationSwitch -Foam::Pstream::floatTransfer -( - "floatTransfer", - 0 -); - +// Allocate a serial communicator. This gets overwritten in parallel mode +// (by Pstream::setParRun()) +Foam::Pstream::communicator serialComm(-1, Foam::labelList(1, 0), false); // Number of processors at which the reduce algorithm changes from linear to // tree @@ -642,7 +913,7 @@ const Foam::debug::optimisationSwitch Foam::Pstream::nProcsSimpleSum ( "nProcsSimpleSum", - 16 + 0 ); @@ -650,11 +921,24 @@ Foam::debug::optimisationSwitch Foam::Pstream::defaultCommsType ( "commsType", -// "nonBlocking", -// "scheduled", - "blocking", + "nonBlocking", "blocking, nonBlocking, scheduled" ); +const Foam::debug::optimisationSwitch +Foam::Pstream::nPollProcInterfaces +( + "nPollProcInterfaces", + 0 +); + +// Default communicator +Foam::label Foam::Pstream::worldComm(0); + + +// Warn for use of any communicator +Foam::label Foam::Pstream::warnComm(-1); + + // ************************************************************************* // diff --git a/src/foam/db/IOstreams/Pstreams/Pstream.H b/src/foam/db/IOstreams/Pstreams/Pstream.H index 59239ed17..2f44eca57 100644 --- a/src/foam/db/IOstreams/Pstreams/Pstream.H +++ b/src/foam/db/IOstreams/Pstreams/Pstream.H @@ -29,11 +29,11 @@ Description SourceFiles Pstream.C - PstreamsPrint.C PstreamCommsStruct.C gatherScatter.C combineGatherScatter.C gatherScatterList.C + PstreamExchange.C \*---------------------------------------------------------------------------*/ @@ -47,6 +47,8 @@ SourceFiles #include "NamedEnum.H" #include "dynamicLabelList.H" #include "optimisationSwitch.H" +#include "ListOps.H" +#include "LIFOStack.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -59,7 +61,6 @@ namespace Foam class Pstream { - public: //- Types of communications @@ -72,6 +73,7 @@ public: static const NamedEnum commsTypeNames; + // Public classes //- Structure for communicating between processors @@ -85,7 +87,7 @@ public: //- procIDs of processors directly below me labelList below_; - //- procIDs of all processors below (so not just directly below) + //- procIDs of all processors below (not just directly below) labelList allBelow_; //- procIDs of all processors not below. (inverse set of @@ -157,7 +159,7 @@ public: friend Ostream& operator<<(Ostream&, const commsStruct&); }; - //- combineReduce operator for lists. Used for counting. + //- combineReduce operator for lists. Used for counting. class listEq { public: @@ -165,7 +167,7 @@ public: template void operator()(T& x, const T& y) const { - forAll(y, i) + forAll (y, i) { if (y[i].size()) { @@ -183,32 +185,38 @@ private: //- Is this a parallel run? static bool parRun_; - //- My processor index - static int myProcNo_; - - //- Process IDs - static List procIDs_; - - //- Default message type + //- Default message type info static const int msgType_; + //- Stack of free comms + static LIFOStack