WIP: code ideas for load/unload functions in dlLibraryTable

.gitignore

@@ -12,7 +12,7 @@
 # File-browser settings - anywhere
 .directory
 
-# Backup/recovery versions - anywhere
+# CVS recovered versions - anywhere
 .#*
 
 # Objects and archives - anywhere
@@ -23,7 +23,7 @@
 # Derived files
 lex.yy.c
 
-# Core dumps
+# Corefiles
 core
 
 # Dependency files - anywhere
@@ -47,21 +47,22 @@ platforms/
 # Reinstate wmake rules that might look like build directories
 !/wmake/rules/*/
 
-# Doxygen generated
+# doxygen generated documentation
 doc/Doxygen/html
 doc/Doxygen/latex
 doc/Doxygen/man
 doc/Doxygen/DTAGS
 
-# Generated files in the main and doc directories
+# Generated files in the main directory (e.g. ReleaseNotes-?.?.html)
+# and in the doc directory
 /*.html
 /doc/*.html
 
-# Untracked configuration/preferences files
+# Untracked configuration files
 /etc/prefs.csh
 /etc/prefs.sh
-/etc/config.csh/prefs.*
-/etc/config.sh/prefs.*
+/etc/config.csh/prefs.csh
+/etc/config.sh/prefs.sh
 /wmake/rules/General/mplibUSER*
 
 # Source packages - anywhere
@@ -72,11 +73,19 @@ doc/Doxygen/DTAGS
 *.tgz
 *.gtgz
 
-# Ignore tags or project files in the main directory
-/.cproject
-/.dir-locals.el
-/.project
+# Ignore the persistent .build tag in the main directory
+/.build
+
+# Ignore .timeStamp in the main directory
+/.timeStamp
+
+# Ignore .tags in the main directory
 /.tags
 
+# Ignore project files in the main directory
+/.cproject
+/.project
+/.dir-locals.el
+
 # Ignore the test directory
 /tutorialsTest

.gitlab/issue_templates/bug.md

@@ -3,7 +3,7 @@
 
   Before opening a new issue, make sure to search for keywords in the issues
   filtered by the "bug" label and check to see if it has already been reported
-
+  
   You can see how your report will be rendered on the platform by using the
   "preview" tab above
 -->

.gitmodules

@@ -4,16 +4,9 @@
 [submodule "avalanche"]
 	path = modules/avalanche
 	url = https://develop.openfoam.com/Community/avalanche.git
+[submodule "catalyst"]
+	path = modules/catalyst
+	url = https://develop.openfoam.com/Community/catalyst.git
 [submodule "adios"]
 	path = modules/adios
 	url = https://develop.openfoam.com/Community/adiosfoam.git
-[submodule "OpenQBMM"]
-	path = modules/OpenQBMM
-	url = https://github.com/OpenQBMM/OpenQBMM.git
-	branch = openfoam.com
-[submodule "visualization"]
-	path = modules/visualization
-	url = https://develop.openfoam.com/modules/visualization.git
-[submodule "external-solver"]
-	path = modules/external-solver
-	url = https://develop.openfoam.com/Modules/external-solver.git

Allwmake

@@ -1,43 +1,38 @@
 #!/bin/sh
 # Run from OPENFOAM top-level directory only
 cd "${0%/*}" || exit
-wmake -check-dir "$WM_PROJECT_DIR" 2>/dev/null || {
+wmakeCheckPwd "$WM_PROJECT_DIR" 2>/dev/null || {
     echo "Error (${0##*/}) : not located in \$WM_PROJECT_DIR"
     echo "    Check your OpenFOAM environment and installation"
     exit 1
 }
-if [ -f "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments ]
-then  . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments || \
-    echo "Argument parse error";
-else
+[ -d "$WM_PROJECT_DIR" -a -f "$WM_PROJECT_DIR/etc/bashrc" ] || {
     echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
     echo "    Check your OpenFOAM environment and installation"
     exit 1
-fi
+}
+
+. "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments
 
 #------------------------------------------------------------------------------
+# Preamble. Report compiler version
+case "$WM_COMPILER" in
+    Gcc*)    gcc --version 2>/dev/null | sed -ne '1p' ;;
+    Clang*)  clang --version 2>/dev/null | sed -ne '1p' ;;
+esac
+
 # Preamble. Report tools or at least the mpirun location
 if [ -f "$WM_PROJECT_DIR"/wmake/scripts/list_tools ]
-then sh "$WM_PROJECT_DIR"/wmake/scripts/list_tools || true
+then  . "$WM_PROJECT_DIR"/wmake/scripts/list_tools ]
 else
     echo "mpirun=$(command -v mpirun || true)"
 fi
-echo
-# Report compiler information
-compiler="$(wmake -show-path-cxx 2>/dev/null || true)"
-if [ -x "$compiler" ]
-then
-    echo "compiler=$compiler"
-    "$compiler" --version 2>/dev/null | sed -ne '1p'
-else
-    echo "compiler=unknown"
-fi
 
 echo
 echo "========================================"
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"
 echo "Starting compile ${WM_PROJECT_DIR##*/} ${0##*/}"
-echo "  $WM_COMPILER ${WM_COMPILER_TYPE:-system} compiler"
+echo "  $WM_COMPILER $WM_COMPILER_TYPE compiler"
 echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
 echo "========================================"
 echo
@@ -48,18 +43,9 @@ echo
 # Compile ThirdParty libraries and applications
 if [ -d "$WM_THIRD_PARTY_DIR" ]
 then
-    if [ -e "$WM_THIRD_PARTY_DIR"/Allwmake.override ]
-    then
-        if [ -x "$WM_THIRD_PARTY_DIR"/Allwmake.override ]
-        then    "$WM_THIRD_PARTY_DIR"/Allwmake.override
-        fi
-    elif [ -x "$WM_THIRD_PARTY_DIR"/Allwmake ]
-    then      "$WM_THIRD_PARTY_DIR"/Allwmake
-    else
-        echo "Skip ThirdParty (no Allwmake* files)"
-    fi
+    "$WM_THIRD_PARTY_DIR/Allwmake"
 else
-    echo "Skip ThirdParty (no directory)"
+    echo "No ThirdParty directory found - skipping"
 fi
 
 echo "========================================"
@@ -73,20 +59,18 @@ echo
 applications/Allwmake $targetType $*
 
 # Additional components/modules
-if [ "$FOAM_MODULE_PREFIX" = false ]
+if [ -d "$WM_PROJECT_DIR/modules" ]
 then
     echo "========================================"
-    echo "OpenFOAM modules disabled (prefix=false)"
+    echo "Compile OpenFOAM modules"
     echo
-elif [ -d "$WM_PROJECT_DIR/modules" ]
-then
     (cd "$WM_PROJECT_DIR/modules" 2>/dev/null && wmake -all)
 fi
 
 # Count files in given directory. Ignore "Test-*" binaries.
 _foamCountDirEntries()
 {
-    (cd "$1" 2>/dev/null && find . -mindepth 1 -maxdepth 1 -type f 2>/dev/null) |\
+    (cd "$1" 2>/dev/null && find -mindepth 1 -maxdepth 1 -type f 2>/dev/null) |\
         sed -e '\@/Test-@d' | wc -l
 }
 
@@ -95,13 +79,13 @@ echo
 date "+%Y-%m-%d %H:%M:%S %z" 2>/dev/null || echo "date is unknown"
 echo "========================================"
 echo "  ${WM_PROJECT_DIR##*/}"
-echo "  $WM_COMPILER ${WM_COMPILER_TYPE:-system} compiler"
+echo "  $WM_COMPILER $WM_COMPILER_TYPE compiler"
 echo "  ${WM_OPTIONS}, with ${WM_MPLIB} ${FOAM_MPI}"
 echo
-echo "  api   = $(etc/openfoam -show-api 2>/dev/null)"
-echo "  patch = $(etc/openfoam -show-patch 2>/dev/null)"
-echo "  bin   = $(_foamCountDirEntries "$FOAM_APPBIN") entries"
-echo "  lib   = $(_foamCountDirEntries "$FOAM_LIBBIN") entries"
+echo "  api   = $(foamEtcFile -show-api 2>/dev/null)"
+echo "  patch = $(foamEtcFile -show-patch 2>/dev/null)"
+echo "  bin   = $(_foamCountDirEntries $FOAM_APPBIN) entries"
+echo "  lib   = $(_foamCountDirEntries $FOAM_LIBBIN) entries"
 echo
 echo "========================================"
 

CONTRIBUTORS.md

@@ -1,56 +0,0 @@
-# Contributors to OpenFOAM
-
-The following is an list of known contributors to OpenFOAM.
-It is likely incomplete...
-
-## Contributors (alphabetical by surname)
-
-- William Bainbridge
-- Gabriel Barajas
-- Kutalmis Bercin
-- Ivor Clifford
-- Greg Collecutt
-- Jonathan Cranford
-- Sergio Ferraris
-- Matej Forman
-- Marian Fuchs
-- Pawan Ghildiyal
-- Chris Greenshields
-- Bernhard Gschaider
-- Andrew Heather
-- David Hill
-- Mattijs Janssens
-- Andrew Jackson
-- Hrvoje Jasak
-- Alexander Kabat vel Job
-- Thilo Knacke
-- Tommaso Lucchini
-- Graham Macpherson
-- Alexey Matveichev
-- Karl Meredith
-- Laurence McGlashan
-- Timo Niemi
-- Haakan Nilsson
-- Niklas Nordin
-- Mark Olesen
-- Victor Olesen
-- Evangelos Papoutsis-Kiachagias
-- Juho Peltola
-- Johan Roenby
-- Henrik Rusche
-- Bruno Santos
-- Henning Scheufler
-- Richard Smith
-- Prashant Sonakar
-- Hilary Spencer
-- Gavin Tabor
-- Zeljko Tukovic
-- Eugene De Villiers
-- Yi Wang
-- Norbert Weber
-- Henry Weller
-- Niklas Wikstrom
-- Thorsten Zirwes
-
-
-<!----------------------------------------------------------------------------->

META-INFO/.gitignore

@@ -5,6 +5,5 @@ build-info
 time-stamp
 
 # Do not track any manifest files
-*[Mm]anifest.txt
-
-####
+Manifest.txt
+manifest.txt

META-INFO/README.md

@@ -1,63 +1,62 @@
-## META-INFO
+# META-INFO
 
-Meta-information is generally for OpenFOAM internal use only.
+Meta-information is for OpenFOAM internal use only.
+
+Do not rely on any files or any file contents in this directory,
+or even the existence of this directory.
 
 The format, content and meaning may be changed at anytime without
-notice. If any of these are changed, these are some of places that
-will need to be updated accordingly:
+notice.
 
-  - bin/foamEtcFile
-  - bin/tools/foamConfigurePaths
-  - bin/tools/foamPackRelease
-  - etc/openfoam
-  - wmake/scripts/wmake-build-info
+The information is provided here for internal documentation purposes.
 
-### api-info
+## api-info
 
 This file and its contents are to be tracked by git.
 
-- File content (api) generated by `wmake -build-info` from the
-  `OPENFOAM` define in `wmake/rules/General/general`
+- File content (api) generated by wmakeBuildInfo from OPENFOAM define
+  in `wmake/rules/General/general`
+
 - File content (patch) is manually generated content.
 
 
-### build-info
+## build-info
 
-This file is ***never*** to be tracked by git, but may be present in
-shipped source archives.
+This file is *never* to be tracked by git, but may be present in shipped
+source archives.
 
-- File content (branch, build) generated by `wmake -build-info` from
-  git information and cached from previous wmake (api)
+- File content (branch, build) generated by wmakeBuildInfo from git
+  information and cached from previous wmake (api)
 
 
-### Content types
+## Content types
 
-#### api
-
-Format: `date +%y%m`
+### api
 
 - 4-digit year-month (YYMM) integer corresponding to the major
   release or in unusual cases an intermediate release.
-  Example, `1712` for the Dec-2017 release.
+
+- Format is year-month, as per `date +%y%m`.
+  Eg, `1712` for the Dec-2017 release.
 
 
-#### patch
-
-Format: `date +%y%m%d`
+### patch
 
 - 6-digit year-month-day (YYMMDD) integer corresponding to a patch-level
   for the given **released** API.
-- The first release can have a patch value of `0` or `1` which
-  indicates that it is unpatched or just released. Alternatively
-  it can have a patch value corresponding to the release date.
+
+- Format is year-month-day, as per `date +%y%m%d`.
+
+- The first release can have a patch value of `0` (unpatched = just
+  released) or a patch value corresponding to the release date.
 
 The patch value is only meaningful together with the api value.
 However, for *development* branches, the patch level should not be
-ascribed too much meaning, but will often correspond to the last
-merge with a *maintenance* (*eg*, `master`) branch.
+ascribed much meaning -- it can be `0` or have a value corresponding
+to the last merge with a *maintenance* (*eg*, `master`) branch.
 
 
-### Flow of information
+## Flow of information
 
 Changes in the build information must be reflected in information
 available in the final binaries. Conversely, it is necessary for later
@@ -71,16 +70,16 @@ distributions to have a record of the same information.
 | build     | git                       | build-info |
 
 
-The command `wmake -build-info -check` is used to determine if
+The command `wmakeBuildInfo -check` is used to determine if
 the saved information needs synchronization. The command
-`wmake -build-info -update` performs the synchronization.
+`wmakeBuildInfo -update` preforms the synchronitzation.
 
 
-### Notes
+## Notes
 
 The saved information is split into two separate files. The `api-info`
 contains more permanent information, whereas the `build-info` is more
 transient in nature.
 
 ----
-2020-06-23
+2019-01-23

META-INFO/api-info

@@ -1,2 +1,2 @@
-api=2012
-patch=210210
+api=1912
+patch=0

README.md

@@ -1,4 +1,4 @@
-## About OpenFOAM
+# About OpenFOAM
 OpenFOAM is a free, open source CFD software [released and developed by OpenCFD Ltd since 2004](http://www.openfoam.com/history/).
 It has a large user base across most areas of engineering and science, from both commercial and academic organisations.
 OpenFOAM has an extensive range of features to solve anything from complex fluid flows involving chemical reactions, turbulence and heat transfer, to acoustics, solid mechanics and electromagnetics.
@@ -11,7 +11,7 @@ individual and group contributors, integrations
 [governance guided activities](https://www.openfoam.com/governance/).
 
 
-## Copyright
+# Copyright
 
 OpenFOAM is free software: you can redistribute it and/or modify it
 under the terms of the GNU General Public License as published by the
@@ -22,7 +22,7 @@ description of the GNU General Public License terms under which you
 may redistribute files.
 
 
-## OpenFOAM Trademark
+# OpenFOAM Trademark
 
 OpenCFD Ltd grants use of its OpenFOAM trademark by Third Parties on a
 licence basis. ESI Group and OpenFOAM Foundation Ltd are currently
@@ -35,7 +35,7 @@ any questions on the use of the OpenFOAM trademark.
 Violations of the Trademark are monitored, and will be duly prosecuted.
 
 
-## Using OpenFOAM
+# Using OpenFOAM
 
 If OpenFOAM has already been compiled on your system, simply source
 the appropriate `etc/bashrc` or `etc/cshrc` file and get started.
@@ -44,7 +44,7 @@ For example, for the OpenFOAM-v1912 version:
 source /installation/path/OpenFOAM-v1912/etc/bashrc
 ```
 
-## Compiling OpenFOAM
+# Compiling OpenFOAM
 
 If you are compiling OpenFOAM from source, please see the relevant
 guides:
@@ -55,7 +55,7 @@ guides:
 | [ThirdParty][repo third] | [readme][link third-readme] | [system requirements][link third-require] | [build][link third-build] |
 
 
-## How do I know which version I am currently using?
+# How do I know which version I am currently using?
 
 The value of the `$WM_PROJECT_DIR` or even `$WM_PROJECT_VERSION` are
 not guaranteed to have any correspondence to the OpenFOAM release
@@ -85,10 +85,10 @@ As can be seen in this example, the git build information is
 supplemented by the date when the last change was authored, which can
 be helpful when the repository contains local changes. If you simply
 wish to know the current API and patch levels directly, the
-`wmake -build-info` provides the relevant information even
+`wmakeBuildInfo` script provides the relevant information even
 when OpenFOAM has not yet been compiled:
 ```
-$ wmake -build-info
+$ wmakeBuildInfo
 make
     api = 1812
     patch = 190531
@@ -110,7 +110,7 @@ patching changes made within the currently active environment and
 should be used with caution.
 
 
-## ThirdParty directory
+# ThirdParty directory
 
 OpenFOAM normally ships with a directory of 3rd-party software and
 build scripts for some 3rd-party software that is either necessary or
@@ -197,14 +197,14 @@ ThirdParty directory will contain either an `Allwmake` file or a
 [link third-require]: https://develop.openfoam.com/Development/ThirdParty-common/blob/develop/Requirements.md
 
 
-## Useful Links
+# Useful Links
 
-- Download [source](https://dl.openfoam.com/source/) and [download and installation instructions](http://www.openfoam.com/download/)
+- Download [source](https://sourceforge.net/projects/openfoam/files/) and [download and installation instructions](http://www.openfoam.com/download/)
 - [Documentation](http://www.openfoam.com/documentation)
 - [Reporting bugs/issues/feature requests](http://www.openfoam.com/code/bug-reporting.php)
-- [Issue tracker](https://develop.openfoam.com/Development/openfoam/-/issues)
-- [Code wiki](https://develop.openfoam.com/Development/openfoam/-/wikis/) and [general wiki](http://wiki.openfoam.com/)
+- [Issue tracker](https://develop.openfoam.com/Development/openfoam/issues)
+- [Wiki](http://wiki.openfoam.com/) and [code wiki](https://develop.openfoam.com/Development/openfoam/wikis/)
 - [Community](http://www.openfoam.com/community/), [Governance](http://www.openfoam.com/governance/)
 - [Contacting OpenCFD](http://www.openfoam.com/contact/)
 
-Copyright 2016-2020 OpenCFD Ltd
+Copyright 2016-2019 OpenCFD Ltd

applications/Allwmake

@@ -1,26 +1,22 @@
 #!/bin/sh
 # Run from OPENFOAM applications/ directory only
 cd "${0%/*}" || exit
-wmake -check-dir "$WM_PROJECT_DIR/applications" 2>/dev/null || {
+wmakeCheckPwd "$WM_PROJECT_DIR/applications" 2>/dev/null || {
     echo "Error (${0##*/}) : not located in \$WM_PROJECT_DIR/applications"
     echo "    Check your OpenFOAM environment and installation"
     exit 1
 }
-if [ -f "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments ]
-then  . "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments || \
-    echo "Argument parse error";
-else
+[ -d "$WM_PROJECT_DIR" -a -f "$WM_PROJECT_DIR/etc/bashrc" ] || {
     echo "Error (${0##*/}) : WM_PROJECT_DIR appears to be incorrect"
     echo "    Check your OpenFOAM environment and installation"
     exit 1
-fi
+}
+
+. "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments
 
 #------------------------------------------------------------------------------
 
 wmake -all $targetType solvers
 wmake -all $targetType utilities
 
-# Optional
-## wmake -all $targetType tools
-
 #------------------------------------------------------------------------------

applications/solvers/DNS/dnsFoam/Allwmake

@@ -2,6 +2,7 @@
 cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/wmake/scripts/AllwmakeParseArguments
 . ${WM_PROJECT_DIR:?}/wmake/scripts/have_fftw
+
 #------------------------------------------------------------------------------
 
 if have_fftw

applications/solvers/DNS/dnsFoam/readTurbulenceProperties.H

@@ -19,3 +19,4 @@
 
     Kmesh K(mesh);
     UOprocess forceGen(K, runTime.deltaTValue(), turbulenceProperties);
+

applications/solvers/acoustic/acousticFoam/Make/files

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

applications/solvers/acoustic/acousticFoam/Make/options

@@ -1,14 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/fvOption/lnInclude \
-    -I$(LIB_SRC)/regionFaModels/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lfvOptions \
-    -lmeshTools \
-    -lsampling \
-    -lregionFaModels

applications/solvers/acoustic/acousticFoam/acousticFoam.C

@@ -1,99 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-   acousticFoam
-
-Group
-    grpAcousticSolvers
-
-Description
-    Acoustic solver solving the acoustic pressure wave equation.
-
-    \f[
-        \ddt2{pa} - c^2 \laplacian{pa} = 0
-    \f]
-
-    where
-    \vartable
-        c       | Sound speed
-        pa      | Acoustic pressure
-    \endvartable
-
-SourceFiles
-    acousticFoam.C
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "fvOptions.H"
-#include "pimpleControl.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    argList::addNote
-    (
-        "Acoustic solver solving the acoustic pressure wave equation."
-    );
-
-    #include "postProcess.H"
-
-    #include "addCheckCaseOptions.H"
-    #include "setRootCaseLists.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
-    #include "createRegionControls.H"
-
-    #include "readTransportProperties.H"
-    #include "createFields.H"
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-        ++runTime;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        while (pimple.correct())
-        {
-            #include "paEqn.H"
-        }
-
-        runTime.write();
-
-        runTime.printExecutionTime(Info);
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/solvers/acoustic/acousticFoam/createFields.H

@@ -1,15 +0,0 @@
-
-Info << "\nReading pa" << endl;
-
-volScalarField pa
-(
-    IOobject
-    (
-        "pa",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);

applications/solvers/acoustic/acousticFoam/createRegionControls.H

@@ -1,8 +0,0 @@
-fvSolution solutionDict(runTime);
-
-const dictionary& pimpleDict = solutionDict.subDict("PIMPLE");
-
-bool solvePrimaryRegion
-(
-    pimpleDict.getOrDefault("solvePrimaryRegion", true)
-);

applications/solvers/acoustic/acousticFoam/paEqn.H

@@ -1,15 +0,0 @@
-
-fvScalarMatrix paEqn
-(
-    fvm::d2dt2(pa) - sqr(c0)*fvc::laplacian(pa)
-);
-
-if (solvePrimaryRegion)
-{
-    paEqn.relax();
-    paEqn.solve();
-}
-else
-{
-    pa.correctBoundaryConditions();
-}

applications/solvers/acoustic/acousticFoam/readTransportProperties.H

@@ -1,23 +0,0 @@
-Info<< "\nReading transportProperties" << endl;
-
-IOdictionary transportProperties
-(
-    IOobject
-    (
-        "transportProperties",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-dimensionedScalar c0("c0", dimVelocity, transportProperties);
-
-dimensionedScalar rho("rho", dimDensity, transportProperties);
-
-scalar MaxCo =
-    max(mesh.surfaceInterpolation::deltaCoeffs()*c0).value()
-   *runTime.deltaT().value();
-
-Info<< "Max acoustic Courant Number = " << MaxCo << endl;

applications/solvers/basic/potentialFoam/createControls.H

@@ -5,5 +5,5 @@ const dictionary& potentialFlow
 
 const int nNonOrthCorr
 (
-    potentialFlow.getOrDefault<int>("nNonOrthogonalCorrectors", 0)
+    potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0)
 );

applications/solvers/basic/potentialFoam/createFields.H

@@ -40,7 +40,7 @@ if (args.found("initialiseUBCs"))
 // converting fixed-value BCs to zero-gradient and vice versa.
 
 // Allow override from command-line -pName option
-const word pName = args.getOrDefault<word>("pName", "p");
+const word pName = args.get<word>("pName", "p");
 
 // Infer the pressure BCs from the velocity
 wordList pBCTypes

applications/solvers/basic/potentialFoam/overPotentialFoam/createControls.H

@@ -5,5 +5,5 @@ const dictionary& potentialFlow
 
 const int nNonOrthCorr
 (
-    potentialFlow.getOrDefault<int>("nNonOrthogonalCorrectors", 0)
+    potentialFlow.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0)
 );

applications/solvers/basic/potentialFoam/overPotentialFoam/createFields.H

@@ -40,7 +40,7 @@ if (args.found("initialiseUBCs"))
 // converting fixed-value BCs to zero-gradient and vice versa.
 
 // Allow override from command-line -pName option
-const word pName = args.getOrDefault<word>("pName", "p");
+const word pName = args.get<word>("pName", "p");
 
 // Infer the pressure BCs from the velocity
 wordList pBCTypes

applications/solvers/combustion/PDRFoam/EaEqn.H

@@ -15,7 +15,7 @@
             )
           : -betav*dpdt
         )
-      - fvm::laplacian(DbThermal, hea)
+      - fvm::laplacian(Db, hea)
       + betav*fvOptions(rho, hea)
     );
 

applications/solvers/combustion/PDRFoam/EauEqn.H

@@ -16,7 +16,13 @@ if (ign.ignited())
             )*rho/thermo.rhou()
           : -betav*dpdt*rho/thermo.rhou()
         )
-      - fvm::laplacian(DbThermal, heau)
+      - fvm::laplacian(Db, heau)
+
+        // These terms cannot be used in partially-premixed combustion due to
+        // the resultant inconsistency between ft and heau transport.
+        // A possible solution would be to solve for ftu as well as ft.
+        //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
+        //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
      ==
         betav*fvOptions(rho, heau)
     );

applications/solvers/combustion/PDRFoam/Make/files

@@ -1,35 +1,33 @@
 XiModels/XiModel/XiModel.C
 XiModels/XiModel/XiModelNew.C
-
-XiModels/transport/transportTwoEqs/transportTwoEqs.C
+XiModels/fixed/fixed.C
+XiModels/algebraic/algebraic.C
+XiModels/transport/transport.C
 
 XiModels/XiEqModels/XiEqModel/XiEqModel.C
 XiModels/XiEqModels/XiEqModel/XiEqModelNew.C
-
 XiModels/XiEqModels/Gulder/Gulder.C
-XiModels/XiEqModels/instability2XiEq/instability2XiEq.C
-XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.C
+XiModels/XiEqModels/instabilityXiEq/instabilityXiEq.C
+XiModels/XiEqModels/SCOPEBlendXiEq/SCOPEBlendXiEq.C
+XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
 
 XiModels/XiGModels/XiGModel/XiGModel.C
 XiModels/XiGModels/XiGModel/XiGModelNew.C
-
 XiModels/XiGModels/KTS/KTS.C
 XiModels/XiGModels/instabilityG/instabilityG.C
-XiModels/XiGModels/instability2G/instability2G.C
 
 PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
 
 PDRModels/dragModels/PDRDragModel/PDRDragModel.C
 PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C
 PDRModels/dragModels/basic/basic.C
-PDRModels/dragModels/basicSch/basicSch.C
 
-PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.C
-
-PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.C
+PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C
+PDRModels/XiGModels/basicXiSubG/basicXiSubG.C
 
 laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
 
+/*PDRFoamAutoRefine.C*/
 PDRFoam.C
 
 EXE = $(FOAM_APPBIN)/PDRFoam

applications/solvers/combustion/PDRFoam/Make/options

@@ -1,9 +1,12 @@
-EXE_INC = -DFULLDEBUG -O0 -g \
+EXE_INC = \
     -IXiModels/XiModel \
     -IXiModels/XiEqModels/XiEqModel \
     -IXiModels/XiGModels/XiGModel \
     -IPDRModels/dragModels/PDRDragModel \
     -IlaminarFlameSpeed/SCOPE \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
@@ -12,13 +15,12 @@ EXE_INC = -DFULLDEBUG -O0 -g \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
     -lmeshTools \
     -lengine \
     -lturbulenceModels \
@@ -28,6 +30,4 @@ EXE_LIBS = \
     -lreactionThermophysicalModels \
     -lspecie \
     -llaminarFlameSpeedModels \
-    -lfiniteVolume \
-    -ldynamicFvMesh \
-    -lfvOptions
+    -ldynamicFvMesh

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,7 +35,7 @@ Description
 
     Combusting RANS code using the b-Xi two-equation model.
     Xi may be obtained by either the solution of the Xi transport
-    equation or from an algebraic expression.  Both approaches are
+    equation or from an algebraic exression.  Both approaches are
     based on Gulder's flame speed correlation which has been shown
     to be appropriate by comparison with the results from the
     spectral model.
@@ -67,7 +66,7 @@ Description
     CR     | Drag tensor (1/m)
     CT     | Turbulence generation parameter (1/m)
     Nv     | Number of obstacles in cell per unit volume (m^-2)
-    nsv    | Tensor whose diagonal indicates the number to subtract from
+    nsv    | Tensor whose diagonal indicates the number to substract from
            | Nv to get the number of obstacles crossing the flow in each
            | direction.
     \endplaintable
@@ -82,11 +81,8 @@ Description
 #include "PDRDragModel.H"
 #include "ignition.H"
 #include "bound.H"
-#include "dynamicFvMesh.H"
-#include "dynamicRefineFvMesh.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -99,24 +95,22 @@ int main(int argc, char *argv[])
     );
 
     #include "postProcess.H"
+
     #include "addCheckCaseOptions.H"
     #include "setRootCaseLists.H"
     #include "createTime.H"
-
-    #include "createDynamicFvMesh.H"
-    #include "createDyMControls.H"
-    #include "initContinuityErrs.H"
-
+    #include "createMesh.H"
+    #include "createControl.H"
     #include "readCombustionProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
-
+    #include "initContinuityErrs.H"
+    #include "createTimeControls.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
     turbulence->validate();
-
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -125,62 +119,33 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readDyMControls.H"
-
+        #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
         ++runTime;
-
         Info<< "\n\nTime = " << runTime.timeName() << endl;
 
+        #include "rhoEqn.H"
+
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            if (pimple.firstIter() || moveMeshOuterCorrectors)
-            {
-                fvc::makeAbsolute(phi, rho, U);
-
-                // Flux estimate for introduced faces.
-                volVectorField rhoU("rhoU", rho*U);
-
-                bool meshChanged = mesh.controlledUpdate();
-
-                if (runTime.write() && meshChanged)
-                {
-                    betav.write();
-                    betai.write();
-                    Blong.write();
-                    Bv.write();
-                    Lobs.write();
-                    CT.write();
-                    drag->writeFields();
-                    flameWrinkling->writeFields();
-                }
-
-                // Make the fluxes relative to the mesh motion
-                fvc::makeRelative(phi, rho, U);
-            }
-
-            if (pimple.firstIter())
-            {
-                #include "rhoEqn.H"
-            }
-
             #include "UEqn.H"
-            #include "bEqn.H"
-            #include "ftEqn.H"
-            #include "EauEqn.H"
-            #include "EaEqn.H"
-
-            if (!ign.ignited())
-            {
-                thermo.heu() == thermo.he();
-            }
 
             // --- Pressure corrector loop
             while (pimple.correct())
             {
+                #include "bEqn.H"
+                #include "ftEqn.H"
+                #include "EauEqn.H"
+                #include "EaEqn.H"
+
+                if (!ign.ignited())
+                {
+                    thermo.heu() == thermo.he();
+                }
+
                 #include "pEqn.H"
             }
 
@@ -190,14 +155,12 @@ int main(int argc, char *argv[])
             }
         }
 
-        rho = thermo.rho();
-
         runTime.write();
 
         runTime.printExecutionTime(Info);
     }
 
-    Info<< "\nEnd\n" << endl;
+    Info<< "\n end\n";
 
     return 0;
 }

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

@@ -0,0 +1,216 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    PDRFoam
+
+Description
+    Solver for compressible premixed/partially-premixed combustion with
+    turbulence modelling.
+
+    Combusting RANS code using the b-Xi two-equation model.
+    Xi may be obtained by either the solution of the Xi transport
+    equation or from an algebraic exression.  Both approaches are
+    based on Gulder's flame speed correlation which has been shown
+    to be appropriate by comparison with the results from the
+    spectral model.
+
+    Strain effects are incorporated directly into the Xi equation
+    but not in the algebraic approximation.  Further work need to be
+    done on this issue, particularly regarding the enhanced removal rate
+    caused by flame compression.  Analysis using results of the spectral
+    model will be required.
+
+    For cases involving very lean Propane flames or other flames which are
+    very strain-sensitive, a transport equation for the laminar flame
+    speed is present.  This equation is derived using heuristic arguments
+    involving the strain time scale and the strain-rate at extinction.
+    the transport velocity is the same as that for the Xi equation.
+
+    For large flames e.g. explosions additional modelling for the flame
+    wrinkling due to surface instabilities may be applied.
+
+    PDR (porosity/distributed resistance) modelling is included to handle
+    regions containing blockages which cannot be resolved by the mesh.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "psiuReactionThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "laminarFlameSpeed.H"
+#include "XiModel.H"
+#include "PDRDragModel.H"
+#include "ignition.H"
+#include "bound.H"
+#include "dynamicRefineFvMesh.H"
+#include "pimpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Solver for compressible premixed/partially-premixed combustion with"
+        " turbulence modelling."
+    );
+
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "readCombustionProperties.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "createTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    turbulence->validate();
+    scalar StCoNum = 0.0;
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    bool hasChanged = false;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        // Indicators for refinement.
+        // Note: before ++runTime only for post-processing reasons.
+        tmp<volScalarField> tmagGradP = mag(fvc::grad(p));
+        volScalarField normalisedGradP
+        (
+            "normalisedGradP",
+            tmagGradP()/max(tmagGradP())
+        );
+        normalisedGradP.writeOpt() = IOobject::AUTO_WRITE;
+        tmagGradP.clear();
+
+        ++runTime;
+
+        Info<< "\n\nTime = " << runTime.timeName() << endl;
+
+        {
+            // Make the fluxes absolute
+            fvc::makeAbsolute(phi, rho, U);
+
+            // Test : disable refinement for some cells
+            bitSet& protectedCell =
+                refCast<dynamicRefineFvMesh>(mesh).protectedCell();
+
+            if (protectedCell.empty())
+            {
+                protectedCell.setSize(mesh.nCells());
+                protectedCell = false;
+            }
+
+            forAll(betav, celli)
+            {
+                if (betav[celli] < 0.99)
+                {
+                    protectedCell.set(celli);
+                }
+            }
+
+            // Flux estimate for introduced faces.
+            volVectorField rhoU("rhoU", rho*U);
+
+            // Do any mesh changes
+            bool meshChanged = mesh.update();
+
+
+            if (meshChanged)
+            {
+                hasChanged = true;
+            }
+
+            if (runTime.write() && hasChanged)
+            {
+                betav.write();
+                Lobs.write();
+                CT.write();
+                drag->writeFields();
+                flameWrinkling->writeFields();
+                hasChanged = false;
+            }
+
+            // Make the fluxes relative to the mesh motion
+            fvc::makeRelative(phi, rho, U);
+        }
+
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "bEqn.H"
+                #include "ftEqn.H"
+                #include "huEqn.H"
+                #include "hEqn.H"
+
+                if (!ign.ignited())
+                {
+                    hu == h;
+                }
+
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "\n end\n";
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicXiSubXiEq.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+    defineTypeNameAndDebug(basicSubGrid, 0);
+    addToRunTimeSelectionTable(XiEqModel, basicSubGrid, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::basicSubGrid::basicSubGrid
+(
+    const dictionary& XiEqProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+
+    B_
+    (
+        IOobject
+        (
+            "B",
+            Su.mesh().facesInstance(),
+            Su.mesh(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        Su.mesh()
+    ),
+
+    XiEqModel_(XiEqModel::New(XiEqModelCoeffs_, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::basicSubGrid::~basicSubGrid()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiEqModels::basicSubGrid::XiEq() const
+{
+    const fvMesh& mesh = Su_.mesh();
+    const volVectorField& U = mesh.lookupObject<volVectorField>("U");
+
+    const volScalarField& Nv = mesh.lookupObject<volScalarField>("Nv");
+    const volSymmTensorField& nsv =
+        mesh.lookupObject<volSymmTensorField>("nsv");
+
+    volScalarField magU(mag(U));
+    volVectorField Uhat
+    (
+        U/(mag(U) + dimensionedScalar("Usmall", U.dimensions(), 1e-4))
+    );
+
+    const scalarField Cw = pow(mesh.V(), 2.0/3.0);
+
+    volScalarField N
+    (
+        IOobject
+        (
+            "N",
+            mesh.time().constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(Nv.dimensions(), Zero)
+    );
+    N.primitiveFieldRef() = Nv.primitiveField()*Cw;
+
+    volSymmTensorField ns
+    (
+        IOobject
+        (
+            "ns",
+            U.mesh().time().timeName(),
+            U.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        U.mesh(),
+        dimensionedSymmTensor(nsv.dimensions(), Zero)
+    );
+    ns.primitiveFieldRef() = nsv.primitiveField()*Cw;
+
+    volScalarField n(max(N - (Uhat & ns & Uhat), scalar(1e-4)));
+    volScalarField b((Uhat & B_ & Uhat)/sqrt(n));
+    volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
+
+    volScalarField XiSubEq
+    (
+        scalar(1)
+      + max(2.2*sqrt(b), min(0.34*magU/up*sqrt(b), scalar(1.6)))
+      * min(n, scalar(1))
+    );
+
+    return (XiSubEq*XiEqModel_->XiEq());
+}
+
+
+bool Foam::XiEqModels::basicSubGrid::read(const dictionary& XiEqProperties)
+{
+    XiEqModel::read(XiEqProperties);
+
+    return XiEqModel_->read(XiEqModelCoeffs_);
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiEqModels::basicSubGrid
+
+Description
+    Basic sub-grid obstacle flame-wrinkling enhancement factor model.
+    Details supplied by J Puttock 2/7/06.
+
+    <b> Sub-grid flame area generation </b>
+
+    \f$ n = N - \hat{\dwea{\vec{U}}}.n_{s}.\hat{\dwea{\vec{U}}} \f$
+    \f$ n_{r} = \sqrt{n} \f$
+
+    where:
+
+        \f$ \hat{\dwea{\vec{U}}} = \dwea{\vec{U}} / \vert \dwea{\vec{U}}
+        \vert \f$
+
+        \f$ b = \hat{\dwea{\vec{U}}}.B.\hat{\dwea{\vec{U}}} / n_{r} \f$
+
+    where:
+
+        \f$ B \f$ is the file "B".
+
+        \f$ N \f$ is the file "N".
+
+        \f$  n_{s} \f$ is the file "ns".
+
+    The flame area enhancement factor \f$ \Xi_{sub} \f$ is expected to
+    approach:
+
+    \f[
+        \Xi_{{sub}_{eq}} =
+            1 + max(2.2 \sqrt{b}, min(0.34 \frac{\vert \dwea{\vec{U}}
+            \vert}{{\vec{U}}^{'}}, 1.6)) \times min(\frac{n}{4}, 1)
+    \f]
+
+
+SourceFiles
+    basicSubGrid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicSubGrid_H
+#define basicSubGrid_H
+
+#include "XiEqModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class basicSubGrid Declaration
+\*---------------------------------------------------------------------------*/
+
+class basicSubGrid
+:
+    public XiEqModel
+{
+    // Private data
+
+        //- tblock
+        volSymmTensorField B_;
+
+        //- Equilibrium Xi model due to turbulence
+        autoPtr<XiEqModel> XiEqModel_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        basicSubGrid(const basicSubGrid&) = delete;
+
+        //- No copy assignment
+        void operator=(const basicSubGrid&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("basicSubGrid");
+
+
+    // Constructors
+
+        //- Construct from components
+        basicSubGrid
+        (
+            const dictionary& XiEqProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su
+        );
+
+
+    //- Destructor
+    virtual ~basicSubGrid();
+
+
+    // Member Functions
+
+        //- Return the flame-wrinkling XiEq
+        virtual tmp<volScalarField> XiEq() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiEqProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiEqModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.C

@@ -1,229 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "normBasicXiSubXiEq.H"
-#include "addToRunTimeSelectionTable.H"
-#include "ignition.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-    defineTypeNameAndDebug(normBasicSubGrid, 0);
-    addToRunTimeSelectionTable(XiEqModel, normBasicSubGrid, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::normBasicSubGrid::normBasicSubGrid
-(
-    const dictionary& XiEqProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiEqModel(XiEqProperties, modelType,thermo, turbulence, Su),
-    Cxpe1_(XiEqModelCoeffs_.get<scalar>("Cxpe1")),
-    Cxpe2_(XiEqModelCoeffs_.get<scalar>("Cxpe2")),
-    Cxpe3_(XiEqModelCoeffs_.get<scalar>("Cxpe3")),
-    Cxpe4_(XiEqModelCoeffs_.get<scalar>("Cxpe4"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::normBasicSubGrid::~normBasicSubGrid()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiEqModels::normBasicSubGrid::XiEq() const
-{
-    const fvMesh& mesh = Su_.mesh();
-    const volVectorField& U = mesh.lookupObject<volVectorField>("U");
-    const volScalarField& b = mesh.lookupObject<volScalarField>("b");
-    const volScalarField& Nv = mesh.lookupObject<volScalarField>("Nv");
-    const volSymmTensorField& nsv =
-        mesh.lookupObject<volSymmTensorField>("nsv");
-    const volSymmTensorField& Bv =
-        mesh.lookupObject<volSymmTensorField>("Bv");
-
-    volScalarField magU(mag(U));
-
-    const scalarField Cw = pow(mesh.V(), 2.0/3.0);
-
-    tmp<volScalarField> tN
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "tN",
-                mesh.time().constant(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar(Nv.dimensions(), Zero)
-        )
-    );
-
-    volScalarField& N = tN.ref();
-
-    N.primitiveFieldRef() = Nv.primitiveField()*Cw;
-
-    tmp<volSymmTensorField> tns
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                "tns",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedSymmTensor(nsv.dimensions(), Zero)
-        )
-    );
-     volSymmTensorField& ns = tns.ref();
-
-    tmp<volSymmTensorField> tB
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                "tB",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedSymmTensor(Bv.dimensions(), Zero)
-        )
-    );
-    volSymmTensorField& B = tB.ref();
-
-    //calculating flame normal
-
-    volVectorField flNormal
-    (
-        "flNormal",
-        fvc::reconstruct(fvc::snGrad(b)*mesh.magSf())
-    );
-
-    volScalarField mgb("mgb", mag(flNormal));
-
-    dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL);
-
-    const volScalarField bc(b*(1.0-b));
-
-    dMgb += 1.0e-8*
-        (bc*mgb)().weightedAverage(mesh.V())
-       /(bc.weightedAverage(mesh.V()) + SMALL);
-
-    mgb += dMgb;
-    flNormal /= mgb;
-
-    B.primitiveFieldRef() = Bv.primitiveField()*sqrt(Cw);
-    volScalarField Ntemp("Ntemp", N);
-    volScalarField Np("Np", max(N - (flNormal & ns & flNormal), scalar(1)));
-
-    // B_ is Bv*sqrt(Cw)
-    volScalarField bl("bl",(flNormal & B & flNormal)/sqrt(Np));
-    bl.min(1.0);
-
-    volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
-
-    IOdictionary combustionProperties
-    (
-        IOobject
-        (
-            "combustionProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-
-    ignition ign(combustionProperties, mesh.time(), mesh);
-
-    dimensionedVector ignLoc("ignLoc", dimLength, ign.sites()[0].location());
-
-    dimensionedScalar filtRad2
-    (
-        "filtRad2",
-        dimLength,
-        6.0*ign.sites()[0].diameter()
-    );
-
-    const volScalarField filDist(mag(mesh.C() - ignLoc));
-
-    const volScalarField filterMult
-    (
-        pos(filDist - filtRad2)*neg(bl - 0.99)*pos(N - 1e-3)
-    );
-
-    tmp<volScalarField> XiSubEq
-    (
-        scalar(1)
-      + min( min(Cxpe1_, Cxpe2_*magU/up)*sqrt(bl), Cxpe3_)
-      * filterMult
-    );
-
-    return XiSubEq;
-}
-
-
-bool Foam::XiEqModels::normBasicSubGrid::read(const dictionary& XiEqProperties)
-{
-    XiEqModel::read(XiEqProperties);
-
-    XiEqModelCoeffs_.readEntry("Cxpe1", Cxpe1_);
-    XiEqModelCoeffs_.readEntry("Cxpe2", Cxpe2_);
-    XiEqModelCoeffs_.readEntry("Cxpe3", Cxpe3_);
-    XiEqModelCoeffs_.readEntry("Cxpe4", Cxpe4_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.H

@@ -1,117 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::XiEqModels::normBasicSubGrid
-
-Description
-
-
-SourceFiles
-    normBasicSubGrid.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef normBasicSubGrid_H
-#define normBasicSubGrid_H
-
-#include "XiEqModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class normBasicSubGrid Declaration
-\*---------------------------------------------------------------------------*/
-
-class normBasicSubGrid
-:
-    public XiEqModel
-{
-    // Private data
-
-        // Constants in the equilibrium Xp equation
-        scalar Cxpe1_;
-        scalar Cxpe2_;
-        scalar Cxpe3_;
-        scalar Cxpe4_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        normBasicSubGrid(const normBasicSubGrid&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const normBasicSubGrid&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("normBasicSubGrid");
-
-
-    // Constructors
-
-        //- Construct from components
-        normBasicSubGrid
-        (
-            const dictionary& XiEqProperties,
-            const word& modelType,
-            const psiuReactionThermo& thermo,
-            const compressible::RASModel& turbulence,
-            const volScalarField& Su
-        );
-
-
-    //- Destructor
-    virtual ~normBasicSubGrid();
-
-
-    // Member Functions
-
-        //- Return the flame-wrinking XiEq
-        virtual tmp<volScalarField> XiEq() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiEqProperties);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace XiEqModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicXiSubG.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiGModels
+{
+    defineTypeNameAndDebug(basicSubGrid, 0);
+    addToRunTimeSelectionTable(XiGModel, basicSubGrid, dictionary);
+};
+};
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiGModels::basicSubGrid::basicSubGrid
+(
+    const dictionary& XiGProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiGModel(XiGProperties, thermo, turbulence, Su),
+    k1(XiGModelCoeffs_.get<scalar>("k1")),
+    XiGModel_(XiGModel::New(XiGModelCoeffs_, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiGModels::basicSubGrid::~basicSubGrid()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiGModels::basicSubGrid::G() const
+{
+    const objectRegistry& db = Su_.db();
+    const volVectorField& U = db.lookupObject<volVectorField>("U");
+    const volScalarField& Nv = db.lookupObject<volScalarField>("Nv");
+    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
+
+    tmp<volScalarField> tGtot = XiGModel_->G();
+    volScalarField& Gtot = tGtot.ref();
+
+    const scalarField Cw = pow(Su_.mesh().V(), 2.0/3.0);
+    scalarField N(Nv.primitiveField()*Cw);
+
+    forAll(N, celli)
+    {
+        if (N[celli] > 1e-3)
+        {
+            Gtot[celli] += k1*mag(U[celli])/Lobs[celli];
+        }
+    }
+
+    return tGtot;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::XiGModels::basicSubGrid::Db() const
+{
+    const objectRegistry& db = Su_.db();
+    const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
+    const volScalarField& rho = db.lookupObject<volScalarField>("rho");
+    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
+    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
+
+    return XiGModel_->Db()
+        + rho*Su_*(Xi - 1.0)*mgb*(0.5*Lobs)*Lobs/(mgb*Lobs + 1.0);
+}
+
+
+bool Foam::XiGModels::basicSubGrid::read(const dictionary& XiGProperties)
+{
+    XiGModel::read(XiGProperties);
+
+    XiGModelCoeffs_.readEntry("k1", k1);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.H

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2013 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiGModels::basicSubGrid
+
+Description
+
+    Basic sub-grid obstacle flame-wrinkling generation rate coefficient model.
+    Details supplied by J Puttock 2/7/06.
+
+    \f$ G_{sub} \f$ denotes the generation coefficient and it is given by
+
+    \f[
+        G_{sub} = k_{1} /frac{\vert \dwea{\vec{U}} \vert}{L_{obs}}
+                 \frac{/Xi_{{sub}_{eq}}-1}{/Xi_{sub}}
+    \f]
+
+    and the removal:
+
+    \f[
+        - k_{1} /frac{\vert \dwea{\vec{U}} \vert}{L_{sub}}
+        \frac{\Xi_{sub}-1}{\Xi_{sub}}
+    \f]
+
+    Finally, \f$ G_{sub} \f$ is added to generation rate \f$ G_{in} \f$
+    due to the turbulence.
+
+SourceFiles
+    basicSubGrid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicSubGrid_H
+#define basicSubGrid_H
+
+#include "XiGModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiGModels
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class basicSubGrid Declaration
+\*---------------------------------------------------------------------------*/
+
+class basicSubGrid
+:
+    public XiGModel
+{
+    // Private data
+
+        //- Sub-grid generation rate coefficient
+        scalar k1;
+
+        //- Xi generation rate model due to turbulence
+        autoPtr<XiGModel> XiGModel_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        basicSubGrid(const basicSubGrid&) = delete;
+
+        //- No copy assignment
+        void operator=(const basicSubGrid&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("basicSubGridG");
+
+
+    // Constructors
+
+        //- Construct from components
+        basicSubGrid
+        (
+            const dictionary& XiGProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su
+        );
+
+
+    //- Destructor
+    virtual ~basicSubGrid();
+
+
+    // Member Functions
+
+        //- Return the flame-wrinkling generation rate
+        virtual tmp<volScalarField> G() const;
+
+        //- Return the flame diffusivity
+        virtual tmp<volScalarField> Db() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiGProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiGModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.C

@@ -1,300 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "normBasicXiSubG.H"
-#include "addToRunTimeSelectionTable.H"
-#include "zeroGradientFvPatchField.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiGModels
-{
-    defineTypeNameAndDebug(normBasicSubGrid, 0);
-    addToRunTimeSelectionTable(XiGModel, normBasicSubGrid, dictionary);
-};
-};
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiGModels::normBasicSubGrid::normBasicSubGrid
-(
-    const dictionary& XiGProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
-//     Bv_
-//     (
-//         IOobject
-//         (
-//             "Bv",
-//             Su.mesh().facesInstance(),
-//             Su.mesh(),
-//             IOobject::MUST_READ,
-//             IOobject::NO_WRITE
-//         ),
-//         Su.mesh()
-//     ),
-    k1_(XiGModelCoeffs_.get<scalar>("k1")),
-    kb1_(XiGModelCoeffs_.get<scalar>("kb1")),
-    kbe_(XiGModelCoeffs_.get<scalar>("kbe")),
-    kbx_(XiGModelCoeffs_.get<scalar>("kbx")),
-    k2_(XiGModelCoeffs_.get<scalar>("k2")),
-    LOverCw_(XiGModelCoeffs_.get<scalar>("LOverCw"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::XiGModels::normBasicSubGrid::~normBasicSubGrid()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::normBasicSubGrid::G() const
-{
-    const objectRegistry& db = Su_.db();
-    const fvMesh& mesh = Su_.mesh();
-    const volVectorField& U = db.lookupObject<volVectorField>("U");
-    const volScalarField& b = db.lookupObject<volScalarField>("b");
-    const volScalarField& Nv = db.lookupObject<volScalarField>("Nv");
-    const volScalarField& St = db.lookupObject<volScalarField>("St");
-    const volSymmTensorField& nsv = db.lookupObject<volSymmTensorField>("nsv");
-    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
-    const volSymmTensorField& Bv = db.lookupObject<volSymmTensorField>("Bv");
-
-    const scalarField Cw(pow(Su_.mesh().V(), 2.0/3.0));
-    volScalarField CwVol
-    (
-        IOobject
-        (
-            "CwVol",
-            mesh.time().timeName(),
-            mesh
-        ),
-        mesh,
-        dimensionSet(dimLength),
-        Cw,
-        zeroGradientFvPatchField<scalar>::typeName
-    );
-    CwVol.correctBoundaryConditions();
-
-    if (!db.foundObject<volScalarField>("Ep"))
-    {
-        FatalErrorIn
-        (
-            "Foam::tmp<Foam::volScalarField> Foam::XiGModels::"
-            "normBasicSubGrid::G() const"
-        )
-        << "Looking for Ep in db which does not exist "
-        << Foam::abort(FatalError);
-    }
-
-    const volScalarField& Ep = db.lookupObject<volScalarField>("Ep");
-    const volScalarField& Xp = db.lookupObject<volScalarField>("Xp");
-    const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
-
-    //tmp<volScalarField> tGtot = XiGModel_->G();
-    tmp<volScalarField> tGtot
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "tGtot",
-                Su_.mesh().time().timeName(),
-                Su_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            Su_.mesh(),
-            dimensionedScalar(inv(dimTime), Zero)
-
-        )
-    );
-
-    volScalarField& Gtot = tGtot.ref();
-
-
-    //calculating flame normal
-
-    volVectorField flNormal(fvc::reconstruct(fvc::snGrad(b)*mesh.magSf()));
-    volScalarField mgb("mgb", mag(flNormal));
-    dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL);
-    {
-        volScalarField bc(b*(1.0-b));
-
-        dMgb += 1.0e-8*
-            (bc*mgb)().weightedAverage(mesh.V())
-           /(bc.weightedAverage(mesh.V()) + SMALL);
-    }
-    mgb += dMgb;
-    flNormal /= mgb;
-
-
-    tmp<volScalarField> tN
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "tN",
-                Su_.mesh().time().timeName(),
-                Su_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            Su_.mesh(),
-            dimensionedScalar(Nv.dimensions(), Zero)
-
-        )
-    );
-
-    volScalarField& N = tN.ref();
-
-    N.primitiveFieldRef() = Nv.primitiveField()*Cw;
-
-    tmp<volSymmTensorField> tns
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                "tns",
-                Su_.mesh().time().timeName(),
-                Su_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            Su_.mesh(),
-            dimensionedSymmTensor(nsv.dimensions(), Zero)
-        )
-    );
-
-    volSymmTensorField& ns = tns.ref();
-
-    ns.primitiveFieldRef() = nsv.primitiveField()*Cw;
-
-    tmp<volSymmTensorField> tB
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                "tB",
-                Su_.mesh().time().timeName(),
-                Su_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            Su_.mesh(),
-            dimensionedSymmTensor(Bv.dimensions(), Zero)
-        )
-    );
-
-    volSymmTensorField& B = tB.ref();
-
-    B.primitiveFieldRef() = Bv.primitiveField()*sqrt(Cw);
-
-    volScalarField Np(max(N - (flNormal & ns & flNormal), scalar(1)));
-
-    // B_ is Bv*sqrt(Cw)
-    volScalarField bl("bl",(flNormal & B & flNormal)/sqrt(Np));
-    bl.min(1.0);
-
-    volScalarField flSpeed("flSpeed", ((U & flNormal) + St)*b/(b+SMALL)) ;
-
-    volScalarField up("up", sqrt((2.0/3.0)*turbulence_.k()));
-
-    const volScalarField Gtot1
-    (
-        "Gtot1",
-        (
-            k1_ + kb1_*min(pow(bl, kbe_), kbx_)
-        )*mag(flSpeed)/(max(Lobs, LOverCw_*CwVol))
-    );
-
-    const volScalarField Gtot2("Gtot2", k2_*Ep*Su_*Xi/(Xp - 0.999));
-
-    const volScalarField value(pos(N - 1.e-3)*neg(bl - 0.99));
-
-    Gtot = value*Gtot1+(1.0 - value)*Gtot2;
-
-
-    //if (Xi.mesh().time().outputTime())
-    {
-        //Gtot.write();
-        //bl.write();
-        //Lobs.write();
-        //flSpeed.write();
-        //N.write();
-    }
-
-    return tGtot;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::normBasicSubGrid::Db() const
-{
-    // Not used //
-    const objectRegistry& db = Su_.db();
-    const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
-    const volScalarField& rho = db.lookupObject<volScalarField>("rho");
-    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
-    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
-    const volScalarField& Db = db.lookupObject<volScalarField>("Db");
-
-    //return  turbulence_.muEff()
-    return  Db
-        + rho*Su_*(Xi - 1.0)*mgb*(0.5*Lobs)*Lobs/(mgb*Lobs + 1.0);
-}
-
-
-bool Foam::XiGModels::normBasicSubGrid::read(const dictionary& XiGProperties)
-{
-    XiGModel::read(XiGProperties);
-
-    XiGModelCoeffs_.readEntry("k1", k1_);
-    XiGModelCoeffs_.readEntry("kb1", kb1_);
-    XiGModelCoeffs_.readEntry("kbe", kbe_);
-    XiGModelCoeffs_.readEntry("kbx", kbx_);
-    XiGModelCoeffs_.readEntry("k2", k2_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.H

@@ -1,132 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::XiGModel::normBasicSubGrid
-
-Description
-
-
-
-SourceFiles
-    normBasicSubGrid.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef normBasicSubGrid_H
-#define normBasicSubGrid_H
-
-#include "XiGModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiGModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class normBasicSubGrid Declaration
-\*---------------------------------------------------------------------------*/
-
-class normBasicSubGrid
-:
-    public XiGModel
-{
-    // Private data
-
-        //- Sub-grid generation rate coefficient
-        scalar k1_;
-
-        //- Sub-grid generation rate coefficient - * sqrt(b)
-        scalar kb1_;
-
-        //- Sub-grid generation rate coefficient - * b
-        scalar kbe_;
-
-        //- Sub-grid generation rate upper limit coefficient - * b
-        scalar kbx_;
-
-        //- Sub-grid generation rate coefficient
-        scalar k2_;
-
-        //- Maximum Lobs/CellWidth
-        scalar LOverCw_;
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        normBasicSubGrid(const normBasicSubGrid&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const normBasicSubGrid&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("normBasicSubGridG");
-
-
-    // Constructors
-
-        //- Construct from components
-        normBasicSubGrid
-        (
-            const dictionary& XiGProperties,
-            const word& modelType,
-            const psiuReactionThermo& thermo,
-            const compressible::RASModel& turbulence,
-            const volScalarField& Su
-        );
-
-
-    //- Destructor
-    virtual ~normBasicSubGrid();
-
-
-    // Member Functions
-
-        //- Return the flame-wrinking generation rate
-        virtual tmp<volScalarField> G() const;
-
-        //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiGProperties);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace XiGModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -71,7 +71,7 @@ Foam::PDRDragModel::PDRDragModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::PDRDragModel::~PDRDragModel()
 {}
@@ -81,7 +81,7 @@ Foam::PDRDragModel::~PDRDragModel()
 
 bool Foam::PDRDragModel::read(const dictionary& PDRProperties)
 {
-    PDRDragModelCoeffs_ = PDRProperties.subDict(type() + "Coeffs");
+    PDRDragModelCoeffs_ = PDRProperties.optionalSubDict(type() + "Coeffs");
 
     PDRDragModelCoeffs_.readEntry("drag", on_);
 

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.H

@@ -76,11 +76,11 @@ private:
 
     // Private Member Functions
 
-        //- Disallow copy construct
-        PDRDragModel(const PDRDragModel&);
+        //- No copy construct
+        PDRDragModel(const PDRDragModel&) = delete;
 
-        //- Disallow default bitwise assignment
-        void operator=(const PDRDragModel&);
+        //- No copy assignment
+        void operator=(const PDRDragModel&) = delete;
 
 
 public:
@@ -161,15 +161,13 @@ public:
         using regIOobject::read;
 
         //- Update properties from given dictionary
-        virtual bool read(const dictionary& PDRProperties);
+        virtual bool read(const dictionary& PDRProperties) = 0;
 
-        //- Write Data
         virtual bool writeData(Ostream&) const
         {
             return true;
         }
 
-        //- Write fields
         virtual void writeFields() const
         {
             NotImplemented;

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,34 +32,32 @@ License
 
 Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New
 (
-    const dictionary& PDRProperties,
+    const dictionary& dict,
     const compressible::RASModel& turbulence,
     const volScalarField& rho,
     const volVectorField& U,
     const surfaceScalarField& phi
 )
 {
-    const word modelType(PDRProperties.get<word>("PDRDragModel"));
+    const word modelType(dict.get<word>("PDRDragModel"));
 
-    Info<< "Selecting flame-wrinkling model " << modelType << endl;
+    Info<< "Selecting drag model " << modelType << endl;
 
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
-        FatalErrorIn
+        FatalIOErrorInLookup
         (
-            "PDRDragModel::New"
-        )   << "Unknown PDRDragModel type "
-            << modelType << nl << nl
-            << "Valid  PDRDragModels are : " << endl
-            << dictionaryConstructorTablePtr_->sortedToc()
-            << exit(FatalError);
+            dict,
+            "PDRDragModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
     }
 
     return autoPtr<PDRDragModel>
-        (cstrIter()(PDRProperties,turbulence, rho, U, phi));
+        (cstrIter()(dict, turbulence, rho, U, phi));
 }
 
 

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2011-2016 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,7 +83,7 @@ Foam::PDRDragModels::basic::basic
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::PDRDragModels::basic::~basic()
 {}
@@ -159,7 +159,6 @@ Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basic::Gk() const
           + Csk*betav*turbulence_.muEff()*sqr(Aw_)*magSqr(U_);
     }
 
-
     return tGk;
 }
 

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.H

@@ -110,11 +110,11 @@ class basic
 
     // Private Member Functions
 
-        //- Disallow copy construct
-        basic(const basic&);
+        //- No copy construct
+        basic(const basic&) = delete;
 
-        //- Disallow default bitwise assignment
-        void operator=(const basic&);
+        //- No copy assignment
+        void operator=(const basic&) = delete;
 
 
 public:

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.C

@@ -1,319 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "basicSch.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace PDRDragModels
-{
-    defineTypeNameAndDebug(basicSch, 0);
-    addToRunTimeSelectionTable(PDRDragModel, basicSch, dictionary);
-}
-}
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::PDRDragModels::basicSch::basicSch
-(
-    const dictionary& PDRProperties,
-    const compressible::RASModel& turbulence,
-    const volScalarField& rho,
-    const volVectorField& U,
-    const surfaceScalarField& phi
-)
-:
-    PDRDragModel(PDRProperties,turbulence, rho, U, phi),
-    Csu("Csu", dimless, PDRDragModelCoeffs_),
-    Csk("Csk", dimless, PDRDragModelCoeffs_),
-    Aw_
-    (
-        IOobject
-        (
-            "Aw",
-            U_.mesh().facesInstance(),
-            U_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh()
-    ),
-
-    CR_
-    (
-        IOobject
-        (
-            "CR",
-            U_.mesh().facesInstance(),
-            U_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh()
-    ),
-    nrCoef_(PDRDragModelCoeffs_.get<scalar>("nrCoef")),
-    nrExp2_(PDRDragModelCoeffs_.get<scalar>("nrExp2")),
-    lCoef_(PDRDragModelCoeffs_.get<scalar>("lCoef")),
-    maxSchFac_(PDRDragModelCoeffs_.get<scalar>("maxSchFac")),
-    subGridSchelkin_(PDRDragModelCoeffs_.get<bool>("subGridSchelkin"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::PDRDragModels::basicSch::~basicSch()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basicSch::Dcu() const
-{
-    tmp<volSymmTensorField> tDragDcu
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                "tDragDcu",
-                U_.mesh().time().constant(),
-                U_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            U_.mesh(),
-            dimensionedSymmTensor(dimMass/dimTime/dimVolume, Zero)
-
-        )
-    );
-
-    volSymmTensorField& DragDcu = tDragDcu.ref();
-
-    if (on_)
-    {
-        const volScalarField& betav =
-            U_.db().lookupObject<volScalarField>("betav");
-
-        DragDcu =
-            (0.5*rho_)*CR_*mag(U_) + (Csu*I)*betav*turbulence_.muEff()*sqr(Aw_);
-    }
-
-    return tDragDcu;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basicSch::Gk() const
-{
-    tmp<volScalarField> tGk
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "tGk",
-                U_.mesh().time().constant(),
-                U_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            U_.mesh(),
-            dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
-
-        )
-    );
-
-    volScalarField& Gk = tGk.ref();
-
-    if (on_)
-    {
-        const volScalarField& betav =
-            U_.db().lookupObject<volScalarField>("betav");
-
-        const volSymmTensorField& CT =
-            U_.db().lookupObject<volSymmTensorField>("CT");
-
-        Gk =
-            (0.5*rho_)*mag(U_)*(U_ & CT & U_)
-          + Csk*betav*turbulence_.muEff()*sqr(Aw_)*magSqr(U_);
-
-        if (subGridSchelkin_)
-        {
-            Gk *= schFac();
-        }
-    }
-
-
-    return tGk;
-}
-
-Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basicSch::schFac() const
-{
-    const volScalarField& Su_ = U_.db().lookupObject<volScalarField>("Su");
-    const volScalarField& rhou_ = U_.db().lookupObject<volScalarField>("rhou");
-    const volScalarField& muu_ = U_.db().lookupObject<volScalarField>("muu");
-
-    tmp<volScalarField> tfac
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "tfac",
-                U_.mesh().time().constant(),
-                U_.mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            U_.mesh(),
-            dimensionedScalar(dimless, Zero)
-
-        )
-    );
-
-    volScalarField& schFac = tfac.ref();
-
-    const volScalarField& k = turbulence_.k();
-    const volScalarField& epsilon = turbulence_.epsilon();
-
-    const volScalarField up(sqrt((2.0/3.0)*k));
-    const volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
-
-    volScalarField Rs(Su_*l*rhou_/muu_);
-
-    if (subGridSchelkin_)
-    {
-        schFac = max
-        (
-            1.0,
-            min
-            (
-                maxSchFac_,
-                pow(Rs, 2.0 * SchelkinExponent(nrCoef_, nrExp2_, Su_))
-            )
-        );
-    }
-
-    return tfac;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basicSch::SchelkinExponent
-(
-    const scalar nrCoef,
-    const scalar nrExp,
-    const volScalarField& Su
-) const
-{
-    const fvMesh& mesh = Su.mesh();
-
-    const volVectorField& U = mesh.lookupObject<volVectorField>("U");
-
-    const volScalarField& Nv = mesh.lookupObject<volScalarField>("Nv");
-    const volSymmTensorField& nsv =
-        mesh.lookupObject<volSymmTensorField>("nsv");
-
-    tmp<volScalarField> tN
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "tN",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedScalar(Nv.dimensions(), Zero)
-
-        )
-    );
-
-    volScalarField& N = tN.ref();
-
-    N.primitiveFieldRef() = Nv.primitiveField()*pow(mesh.V(), 2.0/3.0);
-
-    tmp<volSymmTensorField> tns
-    (
-        new volSymmTensorField
-        (
-            IOobject
-            (
-                "tns",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedSymmTensor(nsv.dimensions(), Zero)
-        )
-    );
-
-    volSymmTensorField& ns = tns.ref();
-
-    ns.primitiveFieldRef() = nsv.primitiveField()*pow(mesh.V(), 2.0/3.0);
-
-    const volVectorField Uhat
-    (
-        U/(mag(U) + dimensionedScalar("Usmall", U.dimensions(), 1e-4))
-    );
-
-    const volScalarField nr(sqrt(max(N - (Uhat & ns & Uhat), scalar(1.0))));
-
-    //Re use tN
-    N.primitiveFieldRef() =
-        nrCoef*((scalar(1.0) - pow(nrExp, nr))/(1.0 - nrExp) - scalar(1.0));
-
-    return tN;
-
-}
-
-
-bool Foam::PDRDragModels::basicSch::read(const dictionary& PDRProperties)
-{
-    PDRDragModel::read(PDRProperties);
-
-    PDRDragModelCoeffs_.readEntry("Csu", Csu.value());
-    PDRDragModelCoeffs_.readEntry("Csk", Csk.value());
-
-    return true;
-}
-
-
-void Foam::PDRDragModels::basicSch::writeFields() const
-{
-    Aw_.write();
-    CR_.write();
-}
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.H

@@ -1,142 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::PDRDragModels::basicSch
-
-Description
-
-
-SourceFiles
-    basicSch.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef basicSch_H
-#define basicSch_H
-
-#include "PDRDragModel.H"
-#include "XiEqModel.H"
-
-
-namespace Foam
-{
-namespace PDRDragModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class basicSch Declaration
-\*---------------------------------------------------------------------------*/
-
-class basicSch
-:
-    public PDRDragModel
-{
-    // Private data
-
-        dimensionedScalar Csu;
-        dimensionedScalar Csk;
-        volScalarField Aw_;
-        volSymmTensorField CR_;
-
-        //- Schelkin effect Model constants
-        const scalar nrCoef_;
-        const scalar nrExp2_;
-        const scalar lCoef_;
-        const scalar maxSchFac_;
-
-        //- Use sub-grid Schelkin effect
-        bool subGridSchelkin_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        basicSch(const basicSch&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const basicSch&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("basicSch");
-
-
-    // Constructors
-
-        //- Construct from components
-        basicSch
-        (
-            const dictionary& PDRProperties,
-            const compressible::RASModel& turbulence,
-            const volScalarField& rho,
-            const volVectorField& U,
-            const surfaceScalarField& phi
-        );
-
-
-    //- Destructor
-    virtual ~basicSch();
-
-
-    // Member Functions
-
-        //- Return the momentum drag coefficient
-        virtual tmp<volSymmTensorField> Dcu() const;
-
-        //- Return the momentum drag turbulence generation rate
-        virtual tmp<volScalarField> Gk() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& PDRProperties);
-
-        //- Write fields
-        void writeFields() const;
-
-        //- Return the schelkin factor for drag turbulence generation rate
-        tmp<volScalarField> schFac() const;
-
-        //- Return the sub-grid Schelkin effect exponent
-        tmp<volScalarField> SchelkinExponent
-        (
-            const scalar,
-            const scalar,
-            const volScalarField&
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace PDRDragModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,7 +29,6 @@ License
 #include "PDRkEpsilon.H"
 #include "PDRDragModel.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Time.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,40 +70,16 @@ PDRkEpsilon::PDRkEpsilon
         modelName
     ),
 
-    C5_(coeffDict_.get<scalar>("C5")),
-    C6_(coeffDict_.get<scalar>("C6")),
-    maxLOverCellW_(coeffDict_.get<scalar>("maxLOverCellW")),
-    lCoef_(coeffDict_.get<scalar>("lCoef")),
-    noTurbUntil_
+    C4_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "noTurbUntil",
+            "C4",
             coeffDict_,
-            0.0
+            0.1
         )
-    ),
-    LOverLobs_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "LOverLobs",
-            coeffDict_,
-            0.2
-        )
-    ),
-    LOverMobs_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "LOverMobs",
-            coeffDict_,
-            0.0
-        )
-    ),
-    bMin_(coeffDict_.get<scalar>("bMin")),
-    flameFilteredG_(coeffDict_.get<bool>("flameFilteredG"))
- {}
+    )
+{}
 
 
 // * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
@@ -119,9 +94,7 @@ bool PDRkEpsilon::read()
 {
     if (RASModel::read())
     {
-        noTurbUntil_.readIfPresent(coeffDict_);
-        LOverLobs_.readIfPresent(coeffDict_);
-        LOverMobs_.readIfPresent(coeffDict_);
+        C4_.readIfPresent(coeffDict_);
         return true;
     }
 
@@ -137,12 +110,16 @@ void PDRkEpsilon::correct()
         nut_ = Cmu_*sqr(k_)/epsilon_;
         nut_.correctBoundaryConditions();
 
+        // Re-calculate thermal diffusivity
+        //***HGWalphat_ = mut_/Prt_;
+        //alphat_.correctBoundaryConditions();
+
         return;
     }
 
     RASModel::correct();
 
-    volScalarField divU("divU",fvc::div(phi_/fvc::interpolate(rho_)));
+    volScalarField divU(fvc::div(phi_/fvc::interpolate(rho_)));
 
     if (mesh_.moving())
     {
@@ -150,9 +127,7 @@ void PDRkEpsilon::correct()
     }
 
     tmp<volTensorField> tgradU = fvc::grad(U_);
-
     volScalarField G(GName(), rho_*nut_*(tgradU() && dev(twoSymm(tgradU()))));
-
     tgradU.clear();
 
     // Update epsilon and G at the wall
@@ -166,102 +141,59 @@ void PDRkEpsilon::correct()
     const volScalarField& Lobs =
         U_.db().lookupObject<volScalarField>("Lobs");
 
-    const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
-
-    const volScalarField& b = mesh_.lookupObject<volScalarField>("b");
-
-    const scalarField Cw(pow(mesh_.V(), 2.0/3.0));
-
     const PDRDragModel& drag =
         U_.db().lookupObject<PDRDragModel>("PDRDragModel");
 
     volScalarField GR(drag.Gk());
 
-    volScalarField LD
+    volScalarField LI
+        (C4_*(Lobs + dimensionedScalar("minLength", dimLength, VSMALL)));
+
+    // Dissipation equation
+    tmp<fvScalarMatrix> epsEqn
     (
-        "LD",
-        (LOverLobs_)*(Lobs + dimensionedScalar("minLength", dimLength, VSMALL))
+        betav*fvm::ddt(rho_, epsilon_)
+      + fvm::div(phi_, epsilon_)
+      - fvm::laplacian(rho_*DepsilonEff(), epsilon_)
+     ==
+        C1_*betav*G*epsilon_/k_
+      + 1.5*pow(Cmu_, 3.0/4.0)*GR*sqrt(k_)/LI
+      - fvm::SuSp(((2.0/3.0)*C1_)*betav*rho_*divU, epsilon_)
+      - fvm::Sp(C2_*betav*rho_*epsilon_/k_, epsilon_)
     );
 
-    LD.primitiveFieldRef() = min(LD.primitiveField(), maxLOverCellW_*sqrt(Cw));
+    epsEqn.ref().relax();
 
-    const volScalarField LI(lCoef_*pow(k_, 3.0/2.0)/epsilon_);
+    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
 
-    IOdictionary combustionProperties
+    solve(epsEqn);
+    bound(epsilon_, epsilonMin_);
+
+
+    // Turbulent kinetic energy equation
+
+    tmp<fvScalarMatrix> kEqn
     (
-        IOobject
-        (
-            "combustionProperties",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
+        betav*fvm::ddt(rho_, k_)
+      + fvm::div(phi_, k_)
+      - fvm::laplacian(rho_*DkEff(), k_)
+     ==
+        betav*G + GR
+      - fvm::SuSp((2.0/3.0)*betav*rho_*divU, k_)
+      - fvm::Sp(betav*rho_*epsilon_/k_, k_)
     );
 
-    ignition ign(combustionProperties, mesh_.time(), U.mesh());
+    kEqn.ref().relax();
+    solve(kEqn);
+    bound(k_, kMin_);
 
-    dimensionedVector ignLoc("ignLoc", dimLength, ign.sites()[0].location());
-
-    dimensionedScalar filtRad2
-    (
-        "filtRad2",
-        dimLength,
-        6.0*ign.sites()[0].diameter()
-    );
-
-    if (flameFilteredG_)
-    {
-        const volScalarField filDist(mag(mesh_.C() - ignLoc));
-        const volScalarField filterG(pos(filDist - filtRad2));
-        const volScalarField filterdivU(pos(filDist - filtRad2));
-        const volScalarField filterGR(pos(filDist - filtRad2)*pos(b - bMin_));
-
-        G *= filterG;
-        GR *= filterGR;
-        divU *= filterdivU;
-    }
-
-    volScalarField Cl(C5_ + (C6_*((LI - LD)/LI)));
-    Cl.max(0.0);
-
-    if (mesh_.time() > noTurbUntil_)
-    {
-        tmp<fvScalarMatrix> epsEqn
-        (
-            betav*fvm::ddt(rho_, epsilon_)
-          + fvm::div(phi_, epsilon_)
-          - fvm::laplacian(rho_*DepsilonEff(), epsilon_)
-        ==
-            C1_*betav*G*epsilon_/k_
-          + Cl*(epsilon_/k_)*GR
-          - fvm::SuSp(((2.0/3.0)*C1_+C3_)*betav*rho_*divU, epsilon_)
-          - fvm::Sp(C2_*betav*rho_*epsilon_/k_, epsilon_)
-        );
-
-        epsEqn.ref().relax();
-        epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
-        solve(epsEqn);
-        bound(epsilon_, epsilonMin_);
-
-        tmp<fvScalarMatrix> kEqn
-        (
-            betav*fvm::ddt(rho_, k_)
-        + fvm::div(phi_, k_)
-        - fvm::laplacian(rho_*DkEff(), k_)
-        ==
-          (betav*G + GR)
-        - fvm::SuSp((2.0/3.0)*betav*rho_*divU, k_)
-        - fvm::Sp(betav*rho_*epsilon_/k_, k_)
-        );
-
-        kEqn.ref().relax();
-        solve(kEqn);
-        bound(k_, kMin_);
-    }
     // Re-calculate viscosity
     nut_ = Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
+
+    // Re-calculate thermal diffusivity
+    //***HGWalphat_ = mut_/Prt_;
+    //alphat_.correctBoundaryConditions();
 }
 
 

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,6 +30,30 @@ Description
     Standard k-epsilon turbulence model with additional source terms
     corresponding to PDR basic drag model (\link basic.H \endlink)
 
+    The default model coefficients correspond to the following:
+    @verbatim
+        PDRkEpsilonCoeffs
+        {
+            Cmu         0.09;
+            C1          1.44;
+            C2          1.92;
+            C3          -0.33;  // only for compressible
+            C4          0.1;
+            sigmak      1.0;    // only for compressible
+            sigmaEps    1.3;
+            Prt         1.0;    // only for compressible
+        }
+    @endverbatim
+
+    The turbulence source term \f$ G_{R} \f$ appears in the
+    \f$ \kappa-\epsilon \f$ equation for the generation of turbulence due to
+    interaction with unresolved obstacles.
+
+    In the \f$ \epsilon  \f$ equation \f$ C_{1} G_{R} \f$ is added as a source
+    term.
+
+    In the \f$ \kappa \f$ equation \f$ G_{R} \f$ is added as a source term.
+
 SourceFiles
     PDRkEpsilon.C
 
@@ -40,7 +63,6 @@ SourceFiles
 #define compressiblePDRkEpsilon_H
 
 #include "kEpsilon.H"
-#include "ignition.H"
 #include "turbulentFluidThermoModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -66,18 +88,10 @@ class PDRkEpsilon
         >
     >
 {
-    // Private Data
+    // Private data
 
         // Model coefficients
-        scalar C5_;
-        scalar C6_;
-        scalar maxLOverCellW_;
-        scalar lCoef_;
-        dimensionedScalar noTurbUntil_;
-        dimensionedScalar LOverLobs_;
-        dimensionedScalar LOverMobs_;
-        scalar bMin_;
-        bool flameFilteredG_;
+        dimensionedScalar C4_;
 
 
 public:

applications/solvers/combustion/PDRFoam/StCourantNo.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2011-2016 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,9 +36,10 @@ Description
 
     if (mesh.nInternalFaces())
     {
-        scalarField sumPhi
+       scalarField sumPhi
         (
-            fvc::surfaceSum(mag(phiSt))().primitiveField()/rho.primitiveField()
+            fvc::surfaceSum(mag(phiSt))().primitiveField()
+          / rho.primitiveField()
         );
 
         StCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

applications/solvers/combustion/PDRFoam/UEqn.H

@@ -1,6 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
         betav*fvm::ddt(rho, U) + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         betav*rho*g

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.C

@@ -1,261 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "BLMgMaXiEq.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-    defineTypeNameAndDebug(BLMgMaXiEq, 0);
-    addToRunTimeSelectionTable(XiEqModel, BLMgMaXiEq, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::BLMgMaXiEq::BLMgMaXiEq
-(
-    const dictionary& XiEqProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su),
-    kaCoef_(XiEqModelCoeffs_.get<scalar>("kaCoef")),
-    lowK0_(XiEqModelCoeffs_.get<scalar>("lowK0")),
-    lowKg_(XiEqModelCoeffs_.get<scalar>("lowKg")),
-    XiEqCoef_(XiEqModelCoeffs_.get<scalar>("XiEqCoef")),
-    alphaCoefP_(XiEqModelCoeffs_.get<scalar>("alphaCoefP")),
-    betaCoefP_(XiEqModelCoeffs_.get<scalar>("betaCoefP")),
-    alphaCoefN_(XiEqModelCoeffs_.get<scalar>("alphaCoefN")),
-    betaCoefN_(XiEqModelCoeffs_.get<scalar>("betaCoefN")),
-    maLim_(XiEqModelCoeffs_.get<scalar>("maLim")),
-    maLim1_(XiEqModelCoeffs_.get<scalar>("maLim1")),
-    quenchCoef_(XiEqModelCoeffs_.get<scalar>("quenchCoef")),
-    quenchExp_(XiEqModelCoeffs_.get<scalar>("quenchExp")),
-    quenchM_(XiEqModelCoeffs_.get<scalar>("quenchM")),
-    quenchRate1_(XiEqModelCoeffs_.get<scalar>("quenchRate1")),
-    quenchRate2_(XiEqModelCoeffs_.get<scalar>("quenchRate2")),
-    lCoef_(XiEqModelCoeffs_.get<scalar>("lCoef")),
-    SuMin_(0.01*Su.average()),
-    uPrimeCoef_(XiEqModelCoeffs_.get<scalar>("uPrimeCoef")),
-    nrExp_(XiEqModelCoeffs_.get<scalar>("nrExp")),
-    subGridSchelkin_(XiEqModelCoeffs_.get<bool>("subGridSchelkin")),
-    MaModel
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "combustionProperties",
-                Su.mesh().time().constant(),
-                Su.mesh(),
-                IOobject::MUST_READ
-            )
-        ),
-        thermo
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::BLMgMaXiEq::~BLMgMaXiEq()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiEqModels::BLMgMaXiEq::XiEq() const
-{
-    const volScalarField& k = turbulence_.k();
-    const volScalarField& epsilon = turbulence_.epsilon();
-    volScalarField up("up", sqrt((2.0/3.0)*k));
-    if (subGridSchelkin_)
-    {
-        up.primitiveFieldRef() +=
-            calculateSchelkinEffect(uPrimeCoef_, nrExp_);
-    }
-
-    volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
-    volScalarField Rl(up*l*thermo_.rhou()/thermo_.muu());
-
-    volScalarField upBySu("upBySu", up/(Su_ + SuMin_));
-    volScalarField K("K", kaCoef_*upBySu*upBySu/sqrt(Rl));
-    volScalarField Ma("Ma", MaModel.Ma());
-
-    volScalarField regime("regime", MaModel.Ma()*scalar(0.0));
-
-    tmp<volScalarField> tXiEq
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "XiEq",
-                epsilon.time().timeName(),
-                epsilon.db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            epsilon.mesh(),
-            dimensionedScalar(dimless, Zero)
-        )
-    );
-
-    const objectRegistry& db = Su_.db();
-
-    const volScalarField& b = db.lookupObject<volScalarField>("b");
-
-    const volScalarField multiP1(0.0*pos(b - 0.99) + 1.0*neg(b - 0.99));
-    const volScalarField multiP2(1.0*pos(b - 0.01) + 0.0*neg(b - 0.01));
-
-    volScalarField& xieq = tXiEq.ref();
-
-    forAll(xieq, celli)
-    {
-        scalar alpha;
-        scalar beta;
-        scalar gulderMa;
-
-        if (Ma[celli]>= 0)
-        {
-            gulderMa =
-                1.0
-              + (0.1402 - 0.007*Ma[celli])
-              * upBySu[celli]*sqrt(upBySu[celli]/K[celli]);
-
-            regime[celli] = multiP1[celli]*multiP2[celli];
-
-        }
-        else
-        {
-            gulderMa =
-                1.0
-              + (0.005*Ma[celli]*Ma[celli]+0.01*Ma[celli] + 0.125)
-              * upBySu[celli]*sqrt(upBySu[celli]/K[celli]);
-
-             regime[celli] = 2*multiP1[celli]*multiP2[celli];
-        }
-
-
-        if (K[celli] < (lowK0_ + lowKg_*Ma[celli]) )
-        {
-            xieq[celli] = gulderMa;
-        }
-        else
-        {
-            if (Ma[celli] >= 0.0)
-            {
-                alpha = alphaCoefP_*(maLim_ - Ma[celli]);
-                beta = betaCoefP_*(maLim_ - Ma[celli]);
-                regime[celli] = 3*multiP1[celli]*multiP2[celli];
-            }
-            else
-            {
-                alpha = alphaCoefN_*(maLim1_ - Ma[celli]) ;
-                beta = betaCoefN_*(maLim_ + Ma[celli]);
-                regime[celli] = 4*multiP1[celli]*multiP2[celli];
-            }
-            xieq[celli] = XiEqCoef_*alpha*pow(K[celli], beta)*upBySu[celli];
-        }
-
-        if (Ma[celli] > -3.0 && Ma[celli] < 11.0)
-        {
-            scalar K0p8 = quenchCoef_*pow( Ma[celli] - quenchM_, quenchExp_);
-            scalar quenchRate = quenchRate1_ + quenchRate2_*Ma[celli];
-            if (K[celli] > (K0p8 - 0.223/quenchRate))
-            {
-                xieq[celli] *= 0.8*exp(-quenchRate*(K[celli] - K0p8));
-                regime[celli] = 5*multiP1[celli]*multiP2[celli];
-            }
-        }
-    }
-
-    forAll(xieq.boundaryField(), patchi)
-    {
-        scalarField& xieqp = xieq.boundaryFieldRef()[patchi];
-        const scalarField& Kp = K.boundaryField()[patchi];
-        const scalarField& Map = Ma.boundaryField()[patchi];
-        const scalarField& upBySup = upBySu.boundaryField()[patchi];
-
-        forAll(xieqp, facei)
-        {
-            scalar alpha;
-            scalar beta;
-
-            if (Map[facei] > 0.0)
-            {
-                alpha = alphaCoefP_*(maLim_ - Map[facei]);
-                beta = betaCoefP_*(maLim_ - Map[facei]);
-            }
-            else
-            {
-                alpha = alphaCoefN_*(maLim_ - Map[facei]);
-                beta = betaCoefN_*(maLim_ + Map[facei]);
-            }
-            xieqp[facei] =
-                XiEqCoef_*alpha*pow(Kp[facei], beta)*upBySup[facei];
-
-            if (Map[facei] > -3.0 && Map[facei] < 11.0)
-            {
-                scalar K0p8 = quenchCoef_*pow(Map[facei] - quenchM_, quenchExp_);
-                scalar quenchRate = quenchRate1_ + quenchRate2_*Ma[facei];
-
-                if (Kp[facei] > (K0p8 - 0.223/quenchRate))
-                {
-                    xieqp[facei] *= 0.8*exp(-quenchRate*(Kp[facei] - K0p8));
-                }
-            }
-            else
-            {
-                Info<<
-                    "Markstein Number out of range for Quench Formulation" << endl;
-            }
-        }
-    }
-
-    return tXiEq;
-}
-
-
-bool Foam::XiEqModels::BLMgMaXiEq::read(const dictionary& XiEqProperties)
-{
-    XiEqModel::read(XiEqProperties);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.H

@@ -1,148 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::XiEqModel::BLMgMaXiEq
-
-Description
-    Model for XiEq based on Bradley, Lawes and Mansour (2011)
-    Cobustion and Falme, 158, 123 correlation
-    with a linear correction function to give a plausible profile for XiEq.
-    See \link SCOPELaminarFlameSpeed.H \endlink for details on the SCOPE laminar
-    flame speed model.
-
-SourceFiles
-    BLMgMaXiEq.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef BLMgMaXiEq_H
-#define BLMgMaXiEq_H
-
-#include "XiEqModel.H"
-#include "SCOPELaminarFlameSpeed.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class BLMgMaXiEq Declaration
-\*---------------------------------------------------------------------------*/
-
-class BLMgMaXiEq
-:
-    public XiEqModel
-{
-    // Private data
-
-        // Model constants
-
-        scalar kaCoef_;
-        scalar lowK0_;
-        scalar lowKg_;
-        scalar XiEqCoef_;
-        scalar alphaCoefP_;
-        scalar betaCoefP_;
-        scalar alphaCoefN_;
-        scalar betaCoefN_;
-        scalar maLim_;
-        scalar maLim1_;
-        scalar quenchCoef_, quenchExp_, quenchM_;
-        scalar quenchRate1_, quenchRate2_;
-        scalar lCoef_;
-
-        //- Minimum Su
-        dimensionedScalar SuMin_;
-
-        //- Schelkin effect Model constants
-        scalar uPrimeCoef_;
-        scalar nrExp_;
-
-        //- Use sub-grid Schelkin effect
-        bool subGridSchelkin_;
-
-        //- The SCOPE laminar flame speed model used to obtain the
-        //  Marstein number.  Note: the laminar flame speed need not be
-        //  obtained form the same model.
-        laminarFlameSpeedModels::SCOPE MaModel;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        BLMgMaXiEq(const BLMgMaXiEq&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const BLMgMaXiEq&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("BLMgMaXiEq");
-
-
-    // Constructors
-
-        //- Construct from components
-        BLMgMaXiEq
-        (
-            const dictionary& XiEqProperties,
-            const word& modelType,
-            const psiuReactionThermo& thermo,
-            const compressible::RASModel& turbulence,
-            const volScalarField& Su
-        );
-
-
-    //- Destructor
-    virtual ~BLMgMaXiEq();
-
-
-    // Member Functions
-
-        //- Return the flame-wrinking XiEq
-        virtual tmp<volScalarField> XiEq() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiEqProperties);
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace XiEqModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2011-2016 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,25 +45,20 @@ namespace XiEqModels
 Foam::XiEqModels::Gulder::Gulder
 (
     const dictionary& XiEqProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 :
-    XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su),
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
     XiEqCoef_(XiEqModelCoeffs_.get<scalar>("XiEqCoef")),
     SuMin_(0.01*Su.average()),
     uPrimeCoef_(XiEqModelCoeffs_.get<scalar>("uPrimeCoef")),
-    nrExp_(XiEqModelCoeffs_.get<scalar>("nrExp")),
-    subGridSchelkin_
-    (
-        XiEqModelCoeffs_.get<bool>("subGridSchelkin")
-    )
+    subGridSchelkin_(XiEqModelCoeffs_.get<bool>("subGridSchelkin"))
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::XiEqModels::Gulder::~Gulder()
 {}
@@ -78,7 +73,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 
     if (subGridSchelkin_)
     {
-        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_, nrExp_);
+        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
     }
 
     volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));
@@ -102,7 +97,6 @@ bool Foam::XiEqModels::Gulder::read(const dictionary& XiEqProperties)
 
     XiEqModelCoeffs_.readEntry("XiEqCoef", XiEqCoef_);
     XiEqModelCoeffs_.readEntry("uPrimeCoef", uPrimeCoef_);
-    XiEqModelCoeffs_.readEntry("nrExp", nrExp_);
     XiEqModelCoeffs_.readEntry("subGridSchelkin", subGridSchelkin_);
 
     return true;

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.H

@@ -65,7 +65,6 @@ class Gulder
 
         //- Schelkin effect Model constant
         scalar uPrimeCoef_;
-        scalar nrExp_;
 
         //- Use sub-grid Schelkin effect
         bool subGridSchelkin_;
@@ -73,11 +72,11 @@ class Gulder
 
     // Private Member Functions
 
-        //- Disallow copy construct
-        Gulder(const Gulder&);
+        //- No copy construct
+        Gulder(const Gulder&) = delete;
 
-        //- Disallow default bitwise assignment
-        void operator=(const Gulder&);
+        //- No copy assignment
+        void operator=(const Gulder&) = delete;
 
 
 public:
@@ -92,7 +91,6 @@ public:
         Gulder
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -105,7 +103,7 @@ public:
 
     // Member Functions
 
-        //- Return the flame-wrinking XiEq
+        //- Return the flame-wrinkling XiEq
         virtual tmp<volScalarField> XiEq() const;
 
         //- Update properties from given dictionary

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEBlendXiEq/SCOPEBlendXiEq.C

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SCOPEBlendXiEq.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+    defineTypeNameAndDebug(SCOPEBlend, 0);
+    addToRunTimeSelectionTable(XiEqModel, SCOPEBlend, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::SCOPEBlend::SCOPEBlend
+(
+    const dictionary& XiEqProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+    XiEqModelL_
+    (
+        XiEqModel::New
+        (
+            XiEqModelCoeffs_.subDict("XiEqModelL"),
+            thermo,
+            turbulence,
+            Su
+        )
+    ),
+    XiEqModelH_
+    (
+        XiEqModel::New
+        (
+            XiEqModelCoeffs_.subDict("XiEqModelH"),
+            thermo,
+            turbulence,
+            Su
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::SCOPEBlend::~SCOPEBlend()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEBlend::XiEq() const
+{
+    return pow
+    (
+        pow4(1.0/XiEqModelL_->XiEq()) + pow4(1.0/XiEqModelH_->XiEq()),
+        -0.25
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEBlendXiEq/SCOPEBlendXiEq.H

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiEqModels::SCOPEBlend
+
+Description
+    Simple SCOPEBlendXiEq model for XiEq based on SCOPEXiEqs correlation
+    with a linear correction function to give a plausible profile for XiEq.
+    See @link SCOPELaminarFlameSpeed.H @endlink for details on the SCOPE
+    laminar flame speed model.
+
+SourceFiles
+    SCOPEBlend.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SCOPEBlend_H
+#define SCOPEBlend_H
+
+#include "XiEqModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class SCOPEBlend Declaration
+\*---------------------------------------------------------------------------*/
+
+class SCOPEBlend
+:
+    public XiEqModel
+{
+    // Private data
+
+        //- Low turbulence intensity equilibrium Xi model
+        autoPtr<XiEqModel> XiEqModelL_;
+
+        //- High turbulence intensity equilibrium Xi model
+        autoPtr<XiEqModel> XiEqModelH_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        SCOPEBlend(const SCOPEBlend&) = delete;
+
+        //- No copy assignment
+        void operator=(const SCOPEBlend&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("SCOPEBlend");
+
+
+    // Constructors
+
+        //- Construct from components
+        SCOPEBlend
+        (
+            const dictionary& XiEqProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su
+        );
+
+
+    //- Destructor
+    virtual ~SCOPEBlend();
+
+
+    // Member Functions
+
+        //- Return the flame-wrinkling XiEq
+        virtual tmp<volScalarField> XiEq() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiEqProperties)
+        {
+            return true;
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiEqModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SCOPEXiEq.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+    defineTypeNameAndDebug(SCOPEXiEq, 0);
+    addToRunTimeSelectionTable(XiEqModel, SCOPEXiEq, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::SCOPEXiEq::SCOPEXiEq
+(
+    const dictionary& XiEqProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+    XiEqCoef_(XiEqModelCoeffs_.get<scalar>("XiEqCoef")),
+    XiEqExp_(XiEqModelCoeffs_.get<scalar>("XiEqExp")),
+    lCoef_(XiEqModelCoeffs_.get<scalar>("lCoef")),
+    SuMin_(0.01*Su.average()),
+    uPrimeCoef_(XiEqModelCoeffs_.get<scalar>("uPrimeCoef")),
+    subGridSchelkin_(XiEqModelCoeffs_.get<bool>("subGridSchelkin")),
+    MaModel
+    (
+        Su.mesh().lookupObject<IOdictionary>("combustionProperties"),
+        thermo
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::SCOPEXiEq::~SCOPEXiEq()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
+{
+    const volScalarField& k = turbulence_.k();
+    const volScalarField& epsilon = turbulence_.epsilon();
+
+    volScalarField up(sqrt((2.0/3.0)*k));
+    if (subGridSchelkin_)
+    {
+        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
+    }
+
+    volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
+    volScalarField Rl(up*l*thermo_.rhou()/thermo_.muu());
+
+    volScalarField upBySu(up/(Su_ + SuMin_));
+    volScalarField K(0.157*upBySu/sqrt(Rl));
+    volScalarField Ma(MaModel.Ma());
+
+    tmp<volScalarField> tXiEq
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "XiEq",
+                epsilon.time().timeName(),
+                epsilon.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            epsilon.mesh(),
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    volScalarField& xieq = tXiEq.ref();
+
+    forAll(xieq, celli)
+    {
+        if (Ma[celli] > 0.01)
+        {
+            xieq[celli] =
+                XiEqCoef_*pow(K[celli]*Ma[celli], -XiEqExp_)*upBySu[celli];
+        }
+    }
+
+    volScalarField::Boundary& xieqBf = xieq.boundaryFieldRef();
+
+    forAll(xieq.boundaryField(), patchi)
+    {
+        scalarField& xieqp = xieqBf[patchi];
+        const scalarField& Kp = K.boundaryField()[patchi];
+        const scalarField& Map = Ma.boundaryField()[patchi];
+        const scalarField& upBySup = upBySu.boundaryField()[patchi];
+
+        forAll(xieqp, facei)
+        {
+            if (Ma[facei] > 0.01)
+            {
+                xieqp[facei] =
+                    XiEqCoef_*pow(Kp[facei]*Map[facei], -XiEqExp_)
+                   *upBySup[facei];
+            }
+        }
+    }
+
+    return tXiEq;
+}
+
+
+bool Foam::XiEqModels::SCOPEXiEq::read(const dictionary& XiEqProperties)
+{
+    XiEqModel::read(XiEqProperties);
+
+    XiEqModelCoeffs_.readEntry("XiEqCoef", XiEqCoef_);
+    XiEqModelCoeffs_.readEntry("XiEqExp", XiEqExp_);
+    XiEqModelCoeffs_.readEntry("lCoef", lCoef_);
+    XiEqModelCoeffs_.readEntry("uPrimeCoef", uPrimeCoef_);
+    XiEqModelCoeffs_.readEntry("subGridSchelkin", subGridSchelkin_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.H

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2013 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiEqModels::SCOPEXiEq
+
+Description
+    Simple SCOPEXiEq model for XiEq based on SCOPEXiEqs correlation
+    with a linear correction function to give a plausible profile for XiEq.
+    See \link SCOPELaminarFlameSpeed.H \endlink for details on the SCOPE laminar
+    flame speed model.
+
+SourceFiles
+    SCOPEXiEq.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SCOPEXiEq_H
+#define SCOPEXiEq_H
+
+#include "XiEqModel.H"
+#include "SCOPELaminarFlameSpeed.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class SCOPEXiEq Declaration
+\*---------------------------------------------------------------------------*/
+
+class SCOPEXiEq
+:
+    public XiEqModel
+{
+    // Private data
+
+        // Model constant
+        scalar XiEqCoef_;
+
+        // Model constant
+        scalar XiEqExp_;
+
+        // Model constant
+        scalar lCoef_;
+
+        //- Minimum Su
+        dimensionedScalar SuMin_;
+
+        //- Schelkin effect Model constant
+        scalar uPrimeCoef_;
+
+        //- Use sub-grid Schelkin effect
+        bool subGridSchelkin_;
+
+        //- The SCOPE laminar flame speed model used to obtain the
+        //  Marstein number.  Note: the laminar flame speed need not be
+        //  obtained form the same model.
+        laminarFlameSpeedModels::SCOPE MaModel;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        SCOPEXiEq(const SCOPEXiEq&) = delete;
+
+        //- No copy assignment
+        void operator=(const SCOPEXiEq&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("SCOPEXiEq");
+
+
+    // Constructors
+
+        //- Construct from components
+        SCOPEXiEq
+        (
+            const dictionary& XiEqProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su
+        );
+
+
+    //- Destructor
+    virtual ~SCOPEXiEq();
+
+
+    // Member Functions
+
+        //- Return the flame-wrinkling XiEq
+        virtual tmp<volScalarField> XiEq() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiEqProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiEqModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,7 +41,6 @@ namespace Foam
 Foam::XiEqModel::XiEqModel
 (
     const dictionary& XiEqProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
@@ -51,7 +50,7 @@ Foam::XiEqModel::XiEqModel
     (
         XiEqProperties.subDict
         (
-            XiEqProperties.get<word>(modelType) + "Coeffs"
+            XiEqProperties.get<word>("XiEqModel") + "Coeffs"
         )
     ),
     thermo_(thermo),
@@ -60,7 +59,7 @@ Foam::XiEqModel::XiEqModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::XiEqModel::~XiEqModel()
 {}
@@ -70,8 +69,7 @@ Foam::XiEqModel::~XiEqModel()
 
 bool Foam::XiEqModel::read(const dictionary& XiEqProperties)
 {
-
-    XiEqModelCoeffs_ = XiEqProperties.subDict(type() + "Coeffs");
+    XiEqModelCoeffs_ = XiEqProperties.optionalSubDict(type() + "Coeffs");
 
     return true;
 }
@@ -79,20 +77,22 @@ bool Foam::XiEqModel::read(const dictionary& XiEqProperties)
 
 void Foam::XiEqModel::writeFields() const
 {
+    //***HGW It is not clear why B is written here
+    if (Su_.mesh().foundObject<volSymmTensorField>("B"))
+    {
+        const volSymmTensorField& B =
+            Su_.mesh().lookupObject<volSymmTensorField>("B");
+        B.write();
+    }
 }
 
 
-
-Foam::tmp<Foam::volScalarField>Foam::XiEqModel::calculateSchelkinEffect
-(
-    const scalar uPrimeCoef,
-    const scalar nrExp
-) const
+Foam::tmp<Foam::volScalarField>
+Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
 {
     const fvMesh& mesh = Su_.mesh();
 
     const volVectorField& U = mesh.lookupObject<volVectorField>("U");
-
     const volSymmTensorField& CT = mesh.lookupObject<volSymmTensorField>("CT");
     const volScalarField& Nv = mesh.lookupObject<volScalarField>("Nv");
     const volSymmTensorField& nsv =
@@ -115,30 +115,22 @@ Foam::tmp<Foam::volScalarField>Foam::XiEqModel::calculateSchelkinEffect
             dimensionedScalar(Nv.dimensions(), Zero)
         )
     );
-
     volScalarField& N = tN.ref();
-
     N.primitiveFieldRef() = Nv.primitiveField()*pow(mesh.V(), 2.0/3.0);
 
-    tmp<volSymmTensorField> tns
+    volSymmTensorField ns
     (
-        new volSymmTensorField
+        IOobject
         (
-            IOobject
-            (
-                "tns",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            "tns",
+            mesh.time().timeName(),
             mesh,
-            dimensionedSymmTensor(nsv.dimensions(), Zero)
-        )
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedSymmTensor(nsv.dimensions(), Zero)
     );
-
-    volSymmTensorField& ns = tns.ref();
-
     ns.primitiveFieldRef() = nsv.primitiveField()*pow(mesh.V(), 2.0/3.0);
 
     const volVectorField Uhat
@@ -152,11 +144,13 @@ Foam::tmp<Foam::volScalarField>Foam::XiEqModel::calculateSchelkinEffect
 
     const scalarField upLocal(uPrimeCoef*sqrt((U & CT & U)*cellWidth));
 
-    //Re use tN
-    N.primitiveFieldRef() = upLocal*(max(scalar(1.0), pow(nr, nrExp)) - 1.0);
+    const scalarField deltaUp(upLocal*(max(scalar(1), pow(nr, 0.5)) - 1.0));
+
+    // Re use tN
+    N.primitiveFieldRef() = upLocal*(max(scalar(1), pow(nr, 0.5)) - 1.0);
 
     return tN;
-
 }
 
+
 // ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.H

@@ -28,6 +28,12 @@ Class
 
 Description
     Base-class for all XiEq models used by the b-XiEq combustion model.
+    The available models are :
+        \link basicXiSubXiEq.H \endlink
+        \link Gulder.H \endlink
+        \link instabilityXiEq.H \endlink
+        \link SCOPEBlendXiEq.H \endlink
+        \link SCOPEXiEq.H \endlink
 
 SourceFiles
     XiEqModel.C
@@ -70,15 +76,16 @@ protected:
         //- Laminar burning velocity
         const volScalarField& Su_;
 
+
 private:
 
     // Private Member Functions
 
-        //- Disallow copy construct
-        XiEqModel(const XiEqModel&);
+        //- No copy construct
+        XiEqModel(const XiEqModel&) = delete;
 
-        //- Disallow default bitwise assignment
-        void operator=(const XiEqModel&);
+        //- No copy assignment
+        void operator=(const XiEqModel&) = delete;
 
 
 public:
@@ -96,14 +103,12 @@ public:
             dictionary,
             (
                 const dictionary& XiEqProperties,
-                const word& modelType,
                 const psiuReactionThermo& thermo,
                 const compressible::RASModel& turbulence,
                 const volScalarField& Su
             ),
             (
                 XiEqProperties,
-                modelType,
                 thermo,
                 turbulence,
                 Su
@@ -117,7 +122,6 @@ public:
         static autoPtr<XiEqModel> New
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -130,7 +134,6 @@ public:
         XiEqModel
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -143,15 +146,14 @@ public:
 
     // Member Functions
 
-        //- Return the flame-wrinking XiEq
-        virtual tmp<volScalarField> XiEq() const = 0;
+        //- Return the flame-wrinkling XiEq
+        virtual tmp<volScalarField> XiEq() const
+        {
+            return turbulence_.muEff();
+        }
 
-        // Only used by sperseded modles
-        tmp<volScalarField> calculateSchelkinEffect
-        (
-            const scalar uPrimeCoef,
-            const scalar nrExp
-        ) const;
+        //- Return the sub-grid Schelkin effect
+        tmp<volScalarField> calculateSchelkinEffect(const scalar) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiEqProperties) = 0;

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModelNew.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,42 +32,30 @@ License
 
 Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
 (
-    const dictionary& propDict,
-    const word& modelName,
+    const dictionary& dict,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 {
-    const word modelType(propDict.get<word>(modelName));
+    const word modelType(dict.get<word>("XiEqModel"));
 
-    Info<< "Selecting flame-wrinkling model Eq " << modelType << endl;
+    Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
-        FatalErrorIn
+        FatalIOErrorInLookup
         (
-            "XiEqModel::New"
-            "("
-            "    const psiuReactionThermo& thermo,"
-            "    const word& modelType,"
-            "    const compressible::RASModel& turbulence,"
-            "    const volScalarField& Su"
-            ")"
-        )   << "Unknown XiEqModel type "
-            << modelType << nl << nl
-            << "Valid XiEqModels are : " << endl
-            << dictionaryConstructorTablePtr_->sortedToc()
-            << exit(FatalError);
+            dict,
+            "XiEqModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
     }
 
-    return autoPtr<XiEqModel>(cstrIter()
-    (
-        propDict, modelName, thermo, turbulence, Su)
-    );
+    return autoPtr<XiEqModel>(cstrIter()(dict, thermo, turbulence, Su));
 }
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.C

@@ -1,141 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "instability2XiEq.H"
-#include "addToRunTimeSelectionTable.H"
-#include "IFstream.H"
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-    defineTypeNameAndDebug(instability2XiEq, 0);
-    addToRunTimeSelectionTable(XiEqModel, instability2XiEq, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::instability2XiEq::instability2XiEq
-(
-    const dictionary& XiEqProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su),
-    saModel_
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "combustionProperties",
-                Su.mesh().time().constant(),
-                Su.mesh(),
-                IOobject::MUST_READ
-            )
-        ),
-        thermo
-    ),
-    CIn_(saModel_.CIn()),
-    defaultCIn_(XiEqModelCoeffs_.get<scalar>("defaultCIn")),
-    XiEqInFade_(XiEqModelCoeffs_.get<scalar>("XiEqInFade")),
-    XiEqModel_
-    (
-        XiEqModel::New(XiEqModelCoeffs_, modelType, thermo, turbulence, Su)
-    )
-{
-    if (CIn_ <=  0.0)
-    {
-        CIn_ = defaultCIn_;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::instability2XiEq::~instability2XiEq()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiEqModels::instability2XiEq::XiEq() const
-{
-    IOdictionary combustionProperties
-    (
-        IOobject
-        (
-            "combustionProperties",
-            Su_.mesh().time().constant(),
-            Su_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    ignition ign(combustionProperties, Su_.mesh().time(), Su_.mesh());
-
-    //const scalar ignTim = ign.sites()[0].tmIgn();
-    scalar curTime = Su_.mesh().time().value();
-    scalar deltaT = Su_.mesh().time().deltaTValue();
-    const scalar ignTim = curTime - deltaT - ign.sites()[0].time();
-
-    volScalarField turbXiEq(XiEqModel_->XiEq());
-
-    volScalarField XiEqIn1("XiEqIn1", 0.0*turbXiEq);
-
-    dimensionedScalar CIn("CIn", dimensionSet(0, -2, 1, 0, 0, 0, 0), CIn_);
-    dimensionedScalar ignTm("ignTm", dimTime, ignTim);
-    XiEqIn1 =  exp(CIn*Su_*Su_*ignTm) - 1.0;
-
-    return
-    (
-        1.0 + sqrt(XiEqInFade_*sqr(XiEqIn1) + sqr(turbXiEq - 1.0))
-    );
-}
-
-
-bool Foam::XiEqModels::instability2XiEq::read(const dictionary& XiEqProperties)
-{
-    XiEqModel::read(XiEqProperties);
-
-    XiEqModelCoeffs_.readEntry("defaultCIn", defaultCIn_);
-    XiEqModelCoeffs_.readEntry("XiEqInFade", XiEqInFade_);
-
-    return XiEqModel_->read(XiEqModelCoeffs_);
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.H

@@ -1,133 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::XiEqModels::instability2XiEq
-
-Description
-
-
-SourceFiles
-    instability2XiEq.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef instability2XiEq_H
-#define instability2XiEq_H
-
-#include "laminarFlameSpeed.H"
-#include "SCOPELaminarFlameSpeed.H"
-#include "ignitionSite.H"
-#include "ignition.H"
-#include "Time.H"
-#include "fvMesh.H"
-#include "XiEqModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class instability2 Declaration
-\*---------------------------------------------------------------------------*/
-
-class instability2XiEq
-:
-    public XiEqModel
-{
-    // Private data
-
-        // -Laminar burning speed
-        laminarFlameSpeedModels::SCOPE saModel_;
-
-        //- GIn (initial instability G)divided by Su^2. Read from fuel file
-        scalar CIn_;
-
-        //- Default CIn if not in fuel file
-        scalar defaultCIn_;
-
-        //- Determines how fast XiEqIn fades out as turbulence comes in
-        scalar XiEqInFade_;
-
-        //- Equilibrium Xi model due to all other effects
-        autoPtr<XiEqModel> XiEqModel_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        instability2XiEq(const instability2XiEq&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const instability2XiEq&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("instability2XiEq");
-
-
-    // Constructors
-
-        //- Construct from components
-        instability2XiEq
-        (
-            const dictionary& XiEqProperties,
-            const word& modelType,
-            const psiuReactionThermo& thermo,
-            const compressible::RASModel& turbulence,
-            const volScalarField& Su
-        );
-
-
-    //- Destructor
-    virtual ~instability2XiEq();
-
-
-    // Member Functions
-
-        //- Return the flame-wrinking XiEq
-        virtual tmp<volScalarField> XiEq() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiEqProperties);
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace XiEqModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instabilityXiEq/instabilityXiEq.C

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "instabilityXiEq.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+    defineTypeNameAndDebug(instability, 0);
+    addToRunTimeSelectionTable(XiEqModel, instability, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::instability::instability
+(
+    const dictionary& XiEqProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+    XiEqIn(XiEqModelCoeffs_.get<scalar>("XiEqIn")),
+    XiEqModel_(XiEqModel::New(XiEqModelCoeffs_, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::instability::~instability()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiEqModels::instability::XiEq() const
+{
+    volScalarField turbXiEq(XiEqModel_->XiEq());
+    return XiEqIn/turbXiEq + turbXiEq;
+}
+
+
+bool Foam::XiEqModels::instability::read(const dictionary& XiEqProperties)
+{
+    XiEqModel::read(XiEqProperties);
+
+    XiEqModelCoeffs_.readEntry("XiEqIn", XiEqIn);
+
+    return XiEqModel_->read(XiEqModelCoeffs_);
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instabilityXiEq/instabilityXiEq.H

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiEqModels::instability
+
+Description
+    This is the equilibrium level of the flame wrinkling generated by
+    instability. It is a constant (default 2.5). It is used in
+    @link XiModel.H @endlink.
+
+SourceFiles
+    instability.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef instability_H
+#define instability_H
+
+#include "XiEqModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class instability Declaration
+\*---------------------------------------------------------------------------*/
+
+class instability
+:
+    public XiEqModel
+{
+    // Private data
+
+        //- Equilibrium Xi due to instability only
+        scalar XiEqIn;
+
+        //- Equilibrium Xi model due to all other effects
+        autoPtr<XiEqModel> XiEqModel_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        instability(const instability&) = delete;
+
+        //- No copy assignment
+        void operator=(const instability&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("instability");
+
+
+    // Constructors
+
+        //- Construct from components
+        instability
+        (
+            const dictionary& XiEqProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su
+        );
+
+
+    //- Destructor
+    virtual ~instability();
+
+
+    // Member Functions
+
+        //- Return the flame-wrinkling XiEq
+        virtual tmp<volScalarField> XiEq() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiEqProperties);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiEqModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C

@@ -45,18 +45,17 @@ namespace XiGModels
 Foam::XiGModels::KTS::KTS
 (
     const dictionary& XiGProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 :
-    XiGModel(XiGProperties,modelType, thermo, turbulence, Su),
+    XiGModel(XiGProperties, thermo, turbulence, Su),
     GEtaCoef_(XiGModelCoeffs_.get<scalar>("GEtaCoef"))
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::XiGModels::KTS::~KTS()
 {}
@@ -75,15 +74,6 @@ Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::G() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::Db() const
-{
-    const objectRegistry& db = Su_.db();
-    const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
-    //return turbulence_.muEff();
-    return Db1;
-}
-
-
 bool Foam::XiGModels::KTS::read(const dictionary& XiGProperties)
 {
     XiGModel::read(XiGProperties);

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2012 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,8 +27,8 @@ Class
     Foam::XiGModels::KTS
 
 Description
-    Simple Kolmogorov time-scale (KTS) model for the flame-wrinkling
-    generation rate.
+    Simple Kolmogorov time-scale (KTS) model for the flame-wrinling generation
+    rate.
 
 SourceFiles
     KTS.C
@@ -62,11 +62,11 @@ class KTS
 
     // Private Member Functions
 
-        //- Disallow copy construct
-        KTS(const KTS&);
+        //- No copy construct
+        KTS(const KTS&) = delete;
 
-        //- Disallow default bitwise assignment
-        void operator=(const KTS&);
+        //- No copy assignment
+        void operator=(const KTS&) = delete;
 
 
 public:
@@ -81,7 +81,6 @@ public:
         KTS
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -94,12 +93,9 @@ public:
 
     // Member Functions
 
-        //- Return the flame-wrinking generation rate
+        //- Return the flame-wrinkling generation rate
         virtual tmp<volScalarField> G() const;
 
-        //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const;
-
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiGProperties);
 };

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModel.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,7 +41,6 @@ namespace Foam
 Foam::XiGModel::XiGModel
 (
     const dictionary& XiGProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
@@ -51,7 +50,7 @@ Foam::XiGModel::XiGModel
     (
         XiGProperties.subDict
         (
-            XiGProperties.get<word>(modelType) + "Coeffs"
+            XiGProperties.get<word>("XiGModel") + "Coeffs"
         )
     ),
     thermo_(thermo),
@@ -60,7 +59,7 @@ Foam::XiGModel::XiGModel
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::XiGModel::~XiGModel()
 {}
@@ -70,7 +69,7 @@ Foam::XiGModel::~XiGModel()
 
 bool Foam::XiGModel::read(const dictionary& XiGProperties)
 {
-    XiGModelCoeffs_ = XiGProperties.subDict(type() + "Coeffs");
+    XiGModelCoeffs_ = XiGProperties.optionalSubDict(type() + "Coeffs");
 
     return true;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModel.H

@@ -28,6 +28,9 @@ Class
 
 Description
     Base-class for all Xi generation models used by the b-Xi combustion model.
+    See Technical Report SH/RE/01R for details on the PDR modelling. For details
+    on the use of XiGModel see \link XiModel.H \endlink. The model available is
+    \link instabilityG.H \endlink
 
 SourceFiles
     XiGModel.C
@@ -69,11 +72,11 @@ private:
 
     // Private Member Functions
 
-        //- Disallow copy construct
-        XiGModel(const XiGModel&);
+        //- No copy construct
+        XiGModel(const XiGModel&) = delete;
 
-        //- Disallow default bitwise assignment
-        void operator=(const XiGModel&);
+        //- No copy assignment
+        void operator=(const XiGModel&) = delete;
 
 
 public:
@@ -91,14 +94,12 @@ public:
             dictionary,
             (
                 const dictionary& XiGProperties,
-                const word& modelType,
                 const psiuReactionThermo& thermo,
                 const compressible::RASModel& turbulence,
                 const volScalarField& Su
             ),
             (
                 XiGProperties,
-                modelType,
                 thermo,
                 turbulence,
                 Su
@@ -112,7 +113,6 @@ public:
         static autoPtr<XiGModel> New
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -125,7 +125,6 @@ public:
         XiGModel
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -138,11 +137,14 @@ public:
 
     // Member Functions
 
-        //- Return the flame-wrinking genration rate
+        //- Return the flame-wrinkling generation rate
         virtual tmp<volScalarField> G() const = 0;
 
         //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const = 0;
+        virtual tmp<volScalarField> Db() const
+        {
+            return turbulence_.muEff();
+        }
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiGProperties) = 0;

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModelNew.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,41 +32,30 @@ License
 
 Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
 (
-    const dictionary& propDict,
-    const word& modelName,
+    const dictionary& dict,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 {
-    const word modelType(propDict.get<word>(modelName));
+    const word modelType(dict.get<word>("XiGModel"));
 
-    Info<< "Selecting flame-wrinkling model for G " << modelType << endl;
+    Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
-        FatalErrorIn
+        FatalIOErrorInLookup
         (
-            "XiGModel::New"
-            "("
-            "    const psiuReactionThermo& thermo,"
-            "    const compressible::RASModel& turbulence,"
-            "    const volScalarField& Su"
-            ")"
-        )   << "Unknown XiGModel type "
-            << modelType << nl << nl
-            << "Valid XiGModels are : " << endl
-            << dictionaryConstructorTablePtr_->sortedToc()
-            << exit(FatalError);
+            dict,
+            "XiGModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
     }
 
-    return autoPtr<XiGModel>(cstrIter()
-    (
-        propDict, modelName, thermo, turbulence, Su)
-    );
+    return autoPtr<XiGModel>(cstrIter()(dict, thermo, turbulence, Su));
 }
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.C

@@ -1,171 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "IFstream.H"
-#include "instability2G.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiGModels
-{
-    defineTypeNameAndDebug(instability2G, 0);
-    addToRunTimeSelectionTable(XiGModel, instability2G, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiGModels::instability2G::instability2G
-(
-    const dictionary& XiGProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
-    saModel_
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "combustionProperties",
-                Su.mesh().time().constant(),
-                Su.mesh(),
-                IOobject::MUST_READ
-            )
-        ),
-        thermo
-    ),
-    CIn_(saModel_.CIn()),
-    defaultCIn_(XiGModelCoeffs_.get<scalar>("defaultCIn")),
-    GInFade_(XiGModelCoeffs_.get<scalar>("GInFade")),
-    GInMult_(XiGModelCoeffs_.get<scalar>("GInMult")),
-    lambdaIn_("lambdaIn", XiGModelCoeffs_),
-    XiGModel_
-    (
-        XiGModel::New(XiGModelCoeffs_,modelType,thermo, turbulence, Su)
-    )
-{
-    if (CIn_ <=  0.0)
-    {
-        CIn_ = defaultCIn_;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::XiGModels::instability2G::~instability2G()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::instability2G::G() const
-{
-        IOdictionary combustionProperties
-        (
-            IOobject
-            (
-                "combustionProperties",
-                Su_.mesh().time().constant(),
-                Su_.mesh(),
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE
-            )
-        );
-
-    ignition ign(combustionProperties, Su_.mesh().time(), Su_.mesh());
-
-    scalar curTime = Su_.mesh().time().value();
-    scalar deltaT = Su_.mesh().time().deltaTValue();
-    const scalar ignTim = curTime - deltaT - ign.sites()[0].time();
-
-    volScalarField turbXiG(XiGModel_->G());
-
-    volScalarField GIn("GIn", 0.0*turbXiG);
-
-
-    forAll (GIn, i)
-    {
-        GIn[i] = CIn_*Su_[i]*Su_[i]*exp(CIn_*Su_[i]*Su_[i]*ignTim)*GInMult_;
-    }
-
-    dimensionedScalar CIn("CIn", dimensionSet(0, -2, 1, 0, 0, 0, 0), CIn_);
-    dimensionedScalar ignTm("ignTm", dimTime, ignTim);
-
-    GIn = CIn*Su_*Su_*exp(CIn*Su_*Su_*ignTm)*GInMult_;
-
-    GIn *=
-        GIn
-        /
-        (
-            GIn
-          + GInFade_*turbXiG
-          + dimensionedScalar("GSmall", inv(dimTime), SMALL)
-        );
-
-    return (GIn + turbXiG);
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::instability2G::Db() const
-{
-    const objectRegistry& db = Su_.db();
-    const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
-    const volScalarField& rho = db.lookupObject<volScalarField>("rho");
-    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
-    const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
-
-    //return  turbulence_.muEff()
-    return Db1
-        + rho*Su_*(Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_);
-}
-
-
-bool Foam::XiGModels::instability2G::read(const dictionary& XiGProperties)
-{
-    XiGModel::read(XiGProperties);
-
-    XiGModelCoeffs_.readEntry("defaultCIn", defaultCIn_);
-    XiGModelCoeffs_.readEntry("GInFade", GInFade_);
-    XiGModelCoeffs_.readEntry("GInMult", GInMult_);
-    XiGModelCoeffs_.readEntry("lambdaIn", lambdaIn_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.H

@@ -1,141 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::XiGModels::instability2G
-
-Description
-
-
-SourceFiles
-    instability2G.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef instability2G_H
-#define instability2G_H
-
-#include "laminarFlameSpeed.H"
-#include "SCOPELaminarFlameSpeed.H"
-#include "XiGModel.H"
-#include "ignitionSite.H"
-#include "ignition.H"
-#include "Time.H"
-#include "fvMesh.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiGModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class instability2G Declaration
-\*---------------------------------------------------------------------------*/
-
-class instability2G
-:
-    public XiGModel
-{
-    // Private data
-
-        laminarFlameSpeedModels::SCOPE saModel_;
-
-        // GIn (inituial instability G)divided by Su^2. Read from fuel file
-        scalar CIn_;
-
-        //- Default CIn if not in fuel file
-        scalar defaultCIn_;
-
-        // Determine how fast GIn fades out as turbulence starts
-        scalar  GInFade_;
-
-        // Set GIn large so that XiEq determines Xi value. Son increase byfactor:
-        scalar  GInMult_;
-
-        //- instability2G length-scale
-        dimensionedScalar lambdaIn_;
-
-        //- Xi generation rate model due to all other processes
-        autoPtr<XiGModel> XiGModel_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        instability2G(const instability2G&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const instability2G&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("instability2G");
-
-
-    // Constructors
-
-        //- Construct from components
-        instability2G
-        (
-            const dictionary& XiGProperties,
-            const word& modelType,
-            const psiuReactionThermo& thermo,
-            const compressible::RASModel& turbulence,
-            const volScalarField& Su
-        );
-
-
-    //- Destructor
-    virtual ~instability2G();
-
-
-    // Member Functions
-
-        //- Return the flame-wrinking generation rate
-        virtual tmp<volScalarField> G() const;
-
-        //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiGProperties);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace XiGModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instabilityG/instabilityG.C

@@ -45,20 +45,19 @@ namespace XiGModels
 Foam::XiGModels::instabilityG::instabilityG
 (
     const dictionary& XiGProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 :
-    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
-    GIn_("GIn", XiGModelCoeffs_),
-    lambdaIn_("lambdaIn", XiGModelCoeffs_),
-    XiGModel_(XiGModel::New(XiGModelCoeffs_,modelType,thermo, turbulence, Su))
+    XiGModel(XiGProperties, thermo, turbulence, Su),
+    GIn_("GIn", dimless/dimTime, XiGModelCoeffs_),
+    lambdaIn_("lambdaIn", dimLength, XiGModelCoeffs_),
+    XiGModel_(XiGModel::New(XiGModelCoeffs_, thermo, turbulence, Su))
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::XiGModels::instabilityG::~instabilityG()
 {}
@@ -77,14 +76,11 @@ Foam::tmp<Foam::volScalarField> Foam::XiGModels::instabilityG::Db() const
 {
     const objectRegistry& db = Su_.db();
     const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
-    const volScalarField& Xp = db.lookupObject<volScalarField>("Xp");
     const volScalarField& rho = db.lookupObject<volScalarField>("rho");
     const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
-    const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
 
-   //return  turbulence_.muEff()
-    return Db1
-        + rho*Su_*(Xp*Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_);
+    return XiGModel_->Db()
+        + rho*Su_*(Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_);
 }
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instabilityG/instabilityG.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2012 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,9 +27,10 @@ Class
     Foam::XiGModels::instabilityG
 
 Description
-    Flame-surface instabilityG flame-wrinking generation rate coefficient model
+    Flame-surface instabilityG flame-wrinkling generation rate coefficient model
     used in \link XiModel.H \endlink.
 
+    See Technical Report SH/RE/01R for details on the PDR modelling.
 
 SourceFiles
     instabilityG.C
@@ -58,7 +59,7 @@ class instabilityG
 {
     // Private data
 
-        //- Flame instabilityG wrinling generation rate coefficient
+        //- Flame instabilityG wrinkling generation rate coefficient
         dimensionedScalar GIn_;
 
         //- InstabilityG length-scale
@@ -70,11 +71,11 @@ class instabilityG
 
     // Private Member Functions
 
-        //- Disallow copy construct
-        instabilityG(const instabilityG&);
+        //- No copy construct
+        instabilityG(const instabilityG&) = delete;
 
-        //- Disallow default bitwise assignment
-        void operator=(const instabilityG&);
+        //- No copy assignment
+        void operator=(const instabilityG&) = delete;
 
 
 public:
@@ -89,7 +90,6 @@ public:
         instabilityG
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -102,7 +102,7 @@ public:
 
     // Member Functions
 
-        //- Return the flame-wrinking generation rate
+        //- Return the flame-wrinkling generation rate
         virtual tmp<volScalarField> G() const;
 
         //- Return the flame diffusivity

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -73,23 +73,11 @@ Foam::XiModel::XiModel
             IOobject::AUTO_WRITE
         ),
         b.mesh()
-    ),
-    Xp_
-    (
-        IOobject
-        (
-            "Xp",
-            b.time().timeName(),
-            b.db(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        b.mesh()
-     )
+    )
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::XiModel::~XiModel()
 {}
@@ -99,7 +87,7 @@ Foam::XiModel::~XiModel()
 
 bool Foam::XiModel::read(const dictionary& XiProperties)
 {
-    XiModelCoeffs_ = XiProperties.subDict(type() + "Coeffs");
+    XiModelCoeffs_ = XiProperties.optionalSubDict(type() + "Coeffs");
 
     return true;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.H

@@ -95,8 +95,8 @@ SourceFiles
 #include "IOdictionary.H"
 #include "psiuReactionThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "fvcDiv.H"
 #include "multivariateSurfaceInterpolationScheme.H"
+#include "fvcDiv.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -124,21 +124,19 @@ protected:
         const volScalarField& b_;
         const surfaceScalarField& phi_;
 
-        //- Flame wrinking field
+        //- Flame wrinkling field
         volScalarField Xi_;
 
-        //- Flame wrinking field for the subgrid obstacles
-        volScalarField Xp_;
 
 private:
 
     // Private Member Functions
 
-        //- Disallow copy construct
-        XiModel(const XiModel&);
+        //- No copy construct
+        XiModel(const XiModel&) = delete;
 
-        //- Disallow default bitwise assignment
-        void operator=(const XiModel&);
+        //- No copy assignment
+        void operator=(const XiModel&) = delete;
 
 
 public:
@@ -211,42 +209,30 @@ public:
 
     // Member Functions
 
-        //- Return the flame-wrinking Xi
+        //- Return the flame-wrinkling Xi
         virtual const volScalarField& Xi() const
         {
             return Xi_;
         }
 
-        //- Return the obstacle-flame-wrinking Xp
-        virtual const volScalarField& Xp() const
-        {
-             return Xp_;
-        }
-
         //- Return the flame diffusivity
         virtual tmp<volScalarField> Db() const
         {
-           const objectRegistry& db = Su_.db();
-           const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
-           //return turbulence_.muEff();
-             return Db1;
+            return turbulence_.muEff();
         }
 
         //- Add Xi to the multivariateSurfaceInterpolationScheme table
         //  if required
-        void addXi
+        virtual void addXi
         (
-            multivariateSurfaceInterpolationScheme<scalar>::fieldTable& fields
+            multivariateSurfaceInterpolationScheme<scalar>::fieldTable&
         )
-        {
-            fields.add(Xi_);
-            fields.add(Xp_);
-        }
+        {}
 
-        //- Correct the flame-wrinking Xi
+        //- Correct the flame-wrinkling Xi
         virtual void correct() = 0;
 
-        //- Correct the flame-wrinking Xi using the given convection scheme
+        //- Correct the flame-wrinkling Xi using the given convection scheme
         virtual void correct(const fv::convectionScheme<scalar>&)
         {
             correct();

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModelNew.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +32,7 @@ License
 
 Foam::autoPtr<Foam::XiModel> Foam::XiModel::New
 (
-    const dictionary& propDict,
+    const dictionary& dict,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su,
@@ -40,27 +41,25 @@ Foam::autoPtr<Foam::XiModel> Foam::XiModel::New
     const surfaceScalarField& phi
 )
 {
-    const word modelType(propDict.get<word>("XiModel"));
+    const word modelType(dict.get<word>("XiModel"));
 
     Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(modelType);
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
-        FatalErrorIn
+        FatalIOErrorInLookup
         (
-            "XiModel::New"
-        )   << "Unknown XiModel type "
-            << modelType << nl << nl
-            << "Valid XiModels are : " << endl
-            << dictionaryConstructorTablePtr_->sortedToc()
-            << exit(FatalError);
+            dict,
+            "XiModel",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
     }
 
     return autoPtr<XiModel>
-        (cstrIter()(propDict, thermo, turbulence, Su, rho, b, phi));
+        (cstrIter()(dict, thermo, turbulence, Su, rho, b, phi));
 }
 
 

applications/solvers/combustion/PDRFoam/XiModels/algebraic/algebraic.C

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "algebraic.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+    defineTypeNameAndDebug(algebraic, 0);
+    addToRunTimeSelectionTable(XiModel, algebraic, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiModels::algebraic::algebraic
+(
+    const dictionary& XiProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su,
+    const volScalarField& rho,
+    const volScalarField& b,
+    const surfaceScalarField& phi
+)
+:
+    XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi),
+    XiShapeCoef(XiModelCoeffs_.get<scalar>("XiShapeCoef")),
+    XiEqModel_(XiEqModel::New(XiProperties, thermo, turbulence, Su)),
+    XiGModel_(XiGModel::New(XiProperties, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiModels::algebraic::~algebraic()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiModels::algebraic::Db() const
+{
+    return XiGModel_->Db();
+}
+
+
+void Foam::XiModels::algebraic::correct()
+{
+    volScalarField XiEqEta(XiEqModel_->XiEq());
+    volScalarField GEta(XiGModel_->G());
+
+    volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999));
+
+    volScalarField XiEqStar(R/(R - GEta));
+
+    Xi_ == 1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b_))*(XiEqStar - 1.0);
+}
+
+
+bool Foam::XiModels::algebraic::read(const dictionary& XiProperties)
+{
+    XiModel::read(XiProperties);
+
+    XiModelCoeffs_.readEntry("XiShapeCoef", XiShapeCoef);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/algebraic/algebraic.H

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiModels::algebraic
+
+Description
+    Simple algebraic model for Xi based on Gulders correlation
+    with a linear correction function to give a plausible profile for Xi.
+    See report TR/HGW/10 for details on the Weller two equations model.
+    See \link XiModel.H \endlink for more details on flame wrinkling modelling.
+
+SourceFiles
+    algebraic.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef algebraic_H
+#define algebraic_H
+
+#include "XiModel.H"
+#include "XiEqModel.H"
+#include "XiGModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class algebraic Declaration
+\*---------------------------------------------------------------------------*/
+
+class algebraic
+:
+    public XiModel
+{
+    // Private data
+
+        scalar XiShapeCoef;
+
+        autoPtr<XiEqModel> XiEqModel_;
+        autoPtr<XiGModel> XiGModel_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        algebraic(const algebraic&) = delete;
+
+        //- No copy assignment
+        void operator=(const algebraic&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("algebraic");
+
+
+    // Constructors
+
+        //- Construct from components
+        algebraic
+        (
+            const dictionary& XiProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su,
+            const volScalarField& rho,
+            const volScalarField& b,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~algebraic();
+
+
+    // Member Functions
+
+        //- Return the flame diffusivity
+        virtual tmp<volScalarField> Db() const;
+
+        //- Correct the flame-wrinkling Xi
+        virtual void correct();
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiProperties);
+
+        //- Write fields of the XiEq model
+        virtual void writeFields()
+        {
+            XiEqModel_().writeFields();
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/fixed/fixed.C

@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixed.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+    defineTypeNameAndDebug(fixed, 0);
+    addToRunTimeSelectionTable(XiModel, fixed, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiModels::fixed::fixed
+(
+    const dictionary& XiProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su,
+    const volScalarField& rho,
+    const volScalarField& b,
+    const surfaceScalarField& phi
+)
+:
+    XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiModels::fixed::~fixed()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::XiModels::fixed::read(const dictionary& XiProperties)
+{
+    return XiModel::read(XiProperties);
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/fixed/fixed.H

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiModels::fixed
+
+Description
+    Fixed value model for Xi. See \link XiModel.H \endlink for more details
+    on flame wrinkling modelling.
+
+SourceFiles
+    fixed.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixed_H
+#define fixed_H
+
+#include "XiModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class fixed Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixed
+:
+    public XiModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        fixed(const fixed&) = delete;
+
+        //- No copy assignment
+        void operator=(const fixed&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fixed");
+
+
+    // Constructors
+
+        //- Construct from components
+        fixed
+        (
+            const dictionary& XiProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su,
+            const volScalarField& rho,
+            const volScalarField& b,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~fixed();
+
+
+    // Member Functions
+
+        //- Correct the flame-wrinkling Xi
+        virtual void correct()
+        {}
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiProperties);
+
+        //- Write fields of the XiEq model
+        virtual void writeFields()
+        {}
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/transport/transport.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "transport.H"
+#include "surfaceInterpolate.H"
+#include "fvmDdt.H"
+#include "fvcLaplacian.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+    defineTypeNameAndDebug(transport, 0);
+    addToRunTimeSelectionTable(XiModel, transport, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiModels::transport::transport
+(
+    const dictionary& XiProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su,
+    const volScalarField& rho,
+    const volScalarField& b,
+    const surfaceScalarField& phi
+)
+:
+    XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi),
+    XiShapeCoef(XiModelCoeffs_.get<scalar>("XiShapeCoef")),
+    XiEqModel_(XiEqModel::New(XiProperties, thermo, turbulence, Su)),
+    XiGModel_(XiGModel::New(XiProperties, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiModels::transport::~transport()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiModels::transport::Db() const
+{
+    return XiGModel_->Db();
+}
+
+
+void Foam::XiModels::transport::correct
+(
+    const fv::convectionScheme<scalar>& mvConvection
+)
+{
+    volScalarField XiEqEta(XiEqModel_->XiEq());
+    volScalarField GEta(XiGModel_->G());
+
+    volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999));
+
+    volScalarField XiEqStar(R/(R - GEta));
+
+    volScalarField XiEq
+    (
+        1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b_))*(XiEqStar - 1.0)
+    );
+
+    volScalarField G(R*(XiEq - 1.0)/XiEq);
+
+    const objectRegistry& db = b_.db();
+    const volScalarField& betav = db.lookupObject<volScalarField>("betav");
+    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
+    const surfaceScalarField& phiSt =
+        db.lookupObject<surfaceScalarField>("phiSt");
+    const volScalarField& Db = db.lookupObject<volScalarField>("Db");
+    const surfaceScalarField& nf = db.lookupObject<surfaceScalarField>("nf");
+
+    surfaceScalarField phiXi
+    (
+        "phiXi",
+        phiSt
+      + (
+          - fvc::interpolate(fvc::laplacian(Db, b_)/mgb)*nf
+          + fvc::interpolate(rho_)*fvc::interpolate(Su_*(1.0/Xi_ - Xi_))*nf
+        )
+    );
+
+    solve
+    (
+        betav*fvm::ddt(rho_, Xi_)
+      + mvConvection.fvmDiv(phi_, Xi_)
+      + fvm::div(phiXi, Xi_)
+      - fvm::Sp(fvc::div(phiXi), Xi_)
+     ==
+        betav*rho_*R
+      - fvm::Sp(betav*rho_*(R - G), Xi_)
+    );
+
+    // Correct boundedness of Xi
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    Xi_.max(1.0);
+    Xi_ = min(Xi_, 2.0*XiEq);
+}
+
+
+bool Foam::XiModels::transport::read(const dictionary& XiProperties)
+{
+    XiModel::read(XiProperties);
+
+    XiModelCoeffs_.readEntry("XiShapeCoef", XiShapeCoef);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/transport/transport.H

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiModels::transport
+
+Description
+    Simple transport model for Xi based on Gulders correlation
+    with a linear correction function to give a plausible profile for Xi.
+    See report TR/HGW/10 for details on the Weller two equations model.
+    See \link XiModel.H \endlink for more details on flame wrinkling modelling.
+
+SourceFiles
+    transport.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef transport_H
+#define transport_H
+
+#include "XiModel.H"
+#include "XiEqModel.H"
+#include "XiGModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class transport Declaration
+\*---------------------------------------------------------------------------*/
+
+class transport
+:
+    public XiModel
+{
+    // Private data
+
+        scalar XiShapeCoef;
+
+        autoPtr<XiEqModel> XiEqModel_;
+        autoPtr<XiGModel> XiGModel_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        transport(const transport&) = delete;
+
+        //- No copy assignment
+        void operator=(const transport&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("transport");
+
+
+    // Constructors
+
+        //- Construct from components
+        transport
+        (
+            const dictionary& XiProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su,
+            const volScalarField& rho,
+            const volScalarField& b,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~transport();
+
+
+    // Member Functions
+
+        //- Return the flame diffusivity
+        virtual tmp<volScalarField> Db() const;
+
+        //- Add Xi to the multivariateSurfaceInterpolationScheme table
+        virtual void addXi
+        (
+            multivariateSurfaceInterpolationScheme<scalar>::fieldTable& fields
+        )
+        {
+            fields.add(Xi_);
+        }
+
+        //- Correct the flame-wrinkling Xi
+        virtual void correct()
+        {
+            NotImplemented;
+        }
+
+        //- Correct the flame-wrinkling Xi using the given convection scheme
+        virtual void correct(const fv::convectionScheme<scalar>& mvConvection);
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiProperties);
+
+        //- Write fields of the XiEq model
+        virtual void writeFields()
+        {
+            XiEqModel_().writeFields();
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.C

@@ -1,253 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "transportTwoEqs.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiModels
-{
-    defineTypeNameAndDebug(transportTwoEqs, 0);
-    addToRunTimeSelectionTable(XiModel, transportTwoEqs, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiModels::transportTwoEqs::transportTwoEqs
-(
-    const dictionary& XiProperties,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su,
-    const volScalarField& rho,
-    const volScalarField& b,
-    const surfaceScalarField& phi
-)
-:
-    XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi),
-    XiShapeCoef_(XiModelCoeffs_.get<scalar>("XiShapeCoef")),
-    CpfiDot_(XiModelCoeffs_.get<scalar>("CpfiDot")),
-    CpfiCross_(XiModelCoeffs_.get<scalar>("CpfiCross")),
-    GEtaExp_(XiModelCoeffs_.get<scalar>("GEtaExp")),
-    LOverCw_(XiModelCoeffs_.get<scalar>("LOverCw")),
-    XiEqModel_
-    (
-        XiEqModel::New(XiProperties, "XiEqModel", thermo, turbulence, Su)
-    ),
-    XiGModel_(XiGModel::New(XiProperties, "XiGModel", thermo, turbulence, Su)),
-    XpEqModel_
-    (
-        XiEqModel::New(XiProperties, "XpEqModel", thermo, turbulence, Su)
-    ),
-    XpGModel_
-    (
-        XiGModel::New(XiProperties, "XpGModel", thermo, turbulence, Su)
-    ),
-    Ep_
-    (
-        IOobject
-        (
-            "Ep",
-            b.time().timeName(),
-            b.db(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        b.mesh()
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
-
-Foam::XiModels::transportTwoEqs::~transportTwoEqs()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-
-Foam::tmp<Foam::volScalarField> Foam::XiModels::transportTwoEqs::Db() const
-{
-    return XiGModel_->Db();
-}
-
-
-void Foam::XiModels::transportTwoEqs::correct
-(
-    const fv::convectionScheme<scalar>& mvConvection
-)
-{
-    const volScalarField XiEqEta(XiEqModel_->XiEq());
-    volScalarField GEta(XiGModel_->G());
-
-    GEta *= max( 1.0, exp( GEtaExp_*(1.0 - (Xi_ - 1.0)/(XiEqEta - 0.999)))) ;
-
-    const volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999));
-
-    const volScalarField XiEqStar(R/(R - GEta));
-
-    const volScalarField XiEq
-    (
-        1.0 + (1.0 + (2*XiShapeCoef_)*(0.5 - b_))*(XiEqStar - 1.0)
-    );
-
-    const volScalarField G(R*(XiEq - 1.0)/XiEq);
-
-
-    const objectRegistry& db = b_.db();
-    const volScalarField& betav = db.lookupObject<volScalarField>("betav");
-    const volScalarField& p = db.lookupObject<volScalarField>("p");
-    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
-    const surfaceScalarField& phiSt =
-        db.lookupObject<surfaceScalarField>("phiSt");
-    const volScalarField& Db = db.lookupObject<volScalarField>("Db");
-    const surfaceScalarField& nf = db.lookupObject<surfaceScalarField>("nf");
-
-    surfaceScalarField phiXi
-    (
-        "phiXi",
-        phiSt
-      + (
-          - fvc::interpolate(fvc::laplacian(Db, b_)/mgb)*nf
-          + fvc::interpolate(rho_)
-          * fvc::interpolate(Su_*(1.0/(Xi_*Xp_) - (Xi_*Xp_)))*nf
-        )
-    );
-
-
-    dimensionedScalar zero
-    (
-        "zero",
-        dimensionSet(2, -6, -2, 0, 0, 0, 0),
-        scalar(0.0)
-    );
-
-    const volScalarField Gpfi
-    (
-          CpfiDot_*sqrt(max(fvc::grad(rho_)&fvc::grad(p), zero ))/rho_*b_*(1.0-b_)
-        + CpfiCross_*sqrt(mag(fvc::grad(rho_)^fvc::grad(p) ))/rho_*b_*(1.0-b_)
-    );
-
-    fvScalarMatrix  XiEqn_
-    (
-        betav*fvm::ddt(rho_, Xi_)
-      + mvConvection.fvmDiv(phi_, Xi_)
-      + fvm::div(phiXi, Xi_)
-      - fvm::Sp(fvc::div(phiXi), Xi_)
-     ==
-        betav*rho_*(R + Gpfi )
-      - fvm::Sp(betav*rho_*(R - G), Xi_)
-    );
-
-    XiEqn_.relax();
-    XiEqn_.solve();
-
-    // Correct boundedness of Xi
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    Xi_.max(1.0);
-    Xi_ = min(Xi_, 2.0*XiEq);
-
-    // Calculation of Xp generated by obstacles
-    volScalarField XpEqEta("XpEqEta",XpEqModel_->XiEq());
-
-    const volScalarField GpEta("GpEta", XpGModel_->G());
-
-    const volScalarField Rp("Rp", GpEta*XpEqEta/(XpEqEta - 0.999));
-
-     const volScalarField XpEq
-    (
-        "XpEq",
-        1.0 + (1.0 + (2*XiShapeCoef_)*(0.5 - b_))*(XpEqEta - 1.0)
-        );
-
-    const volScalarField Gpp("Gpp", Rp*(XpEq - 1.0)/XpEq);
-
-
-    fvScalarMatrix XpEqn_
-    (
-        betav*fvm::ddt(rho_, Xp_)
-      + mvConvection.fvmDiv(phi_, Xp_)
-      + fvm::div(phiXi, Xp_)
-      - fvm::Sp(fvc::div(phiXi), Xp_)
-     ==
-        betav*rho_*Rp
-      - fvm::Sp(betav*rho_*(Rp - Gpp), Xp_)
-    );
-
-    XpEqn_.relax();
-    XpEqn_.solve();
-
-    Xp_.max(1.0);
-    Xp_ = min(Xp_, 20.0*XpEq);
-
-    // Calculate Ep
-    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
-    const scalarField Cw = pow(Su_.mesh().V(), 2.0/3.0);
-    volScalarField LI(Lobs);
-
-    LI.primitiveFieldRef() = max(LI.primitiveField(),LOverCw_*sqrt(Cw));
-
-    fvScalarMatrix EpEqn_
-    (
-        betav*fvm::ddt(rho_, Ep_)
-      + mvConvection.fvmDiv(phi_, Ep_)
-      + fvm::div(phiXi, Ep_)
-      - fvm::Sp(fvc::div(phiXi), Ep_)
-     ==
-        betav*rho_*Gpp*Xp_/LI
-      - fvm::Sp(betav*rho_*Rp, Ep_)
-   );
-
-    EpEqn_.relax();
-    EpEqn_.solve();
-
-    Ep_.max(0.0);
-    Ep_.min(100000.0);
-}
-
-
-bool Foam::XiModels::transportTwoEqs::read(const dictionary& XiProperties)
-{
-    XiModel::read(XiProperties);
-
-    XiModelCoeffs_.readEntry("XiShapeCoef", XiShapeCoef_);
-    XiModelCoeffs_.readEntry("CpfiDot", CpfiDot_);
-    XiModelCoeffs_.readEntry("CpfiCross", CpfiCross_);
-    XiModelCoeffs_.readEntry("GEtaExp", GEtaExp_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.H

@@ -1,154 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::XiModels::transportTwoEqs
-
-Description
-
-SourceFiles
-    transportTwoEqs.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef transportTwoEqs_H
-#define transportTwoEqs_H
-
-#include "XiModel.H"
-#include "XiEqModel.H"
-#include "XiGModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class transportTwoEqs Declaration
-\*---------------------------------------------------------------------------*/
-
-class transportTwoEqs
-:
-    public XiModel
-{
-    // Private data
-
-        scalar XiShapeCoef_;
-        scalar CpfiDot_;
-        scalar CpfiCross_;
-        scalar GEtaExp_;
-
-        //- Maximum Lobs/CellWidth
-        scalar LOverCw_;
-
-        //- Equilibrium for Xi (turbulence)
-        autoPtr<XiEqModel> XiEqModel_;
-
-        //- Generation for Xi (turbulence)
-        autoPtr<XiGModel> XiGModel_;
-
-        //- Equilibrium for Xp (obstacles)
-        autoPtr<XiEqModel> XpEqModel_;
-
-        //- Generation for Xp (obstacles)
-        autoPtr<XiGModel> XpGModel_;
-
-        //- Dissipation lenght scale for subgrid obstacles
-        volScalarField Ep_;
-
-
-    // Private Member Functions
-
-        //- Disallow copy construct
-        transportTwoEqs(const transportTwoEqs&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const transportTwoEqs&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("transportTwoEqs");
-
-
-    // Constructors
-
-        //- Construct from components
-        transportTwoEqs
-        (
-            const dictionary& XiProperties,
-            const psiuReactionThermo& thermo,
-            const compressible::RASModel& turbulence,
-            const volScalarField& Su,
-            const volScalarField& rho,
-            const volScalarField& b,
-            const surfaceScalarField& phi
-        );
-
-
-    //- Destructor
-    virtual ~transportTwoEqs();
-
-
-    // Member Functions
-
-        //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const;
-
-        //- Correct the flame-wrinking Xi
-        virtual void correct()
-        {
-            NotImplemented;
-        }
-
-        //- Correct the flame-wrinking Xi using the given convection scheme
-        virtual void correct(const fv::convectionScheme<scalar>& mvConvection);
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiProperties);
-
-        //- Write fields of the XiEq model
-        virtual void writeFields()
-        {
-            XiEqModel_().writeFields();
-            XpEqModel_().writeFields();
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace XiModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/bEqn.H

@@ -5,20 +5,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,ft_b_h_hu)")
+        mesh.divScheme("div(phi,ft_b_ha_hau)")
     )
 );
 
-rhou = thermo.rhou();
-muu = thermo.muu();
-
-volScalarField Db("Db", turbulence->mut());
-const volScalarField DbThermal
-(
-    "DbAlpha",
-    turbulence->alphaEff()
-    //thermo.alpha() + turbulence->mut()/Prt
-);
+volScalarField Db("Db", turbulence->muEff());
 
 if (ign.ignited())
 {
@@ -26,30 +17,30 @@ if (ign.ignited())
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     Su = unstrainedLaminarFlameSpeed()();
 
-    // Wrinkling due to turbulence
     const volScalarField& Xi = flameWrinkling->Xi();
-    // Wrinkling due to obstacles
-    const volScalarField& Xp = flameWrinkling->Xp();
 
     // progress variable
     // ~~~~~~~~~~~~~~~~~
-    c = (1.0 - b);
+    volScalarField c("c", 1.0 - b);
 
+    // Unburnt gas density
+    // ~~~~~~~~~~~~~~~~~~~
+    volScalarField rhou(thermo.rhou());
 
     // Calculate flame normal etc.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+    //volVectorField n(fvc::grad(b));
     volVectorField n(fvc::reconstruct(fvc::snGrad(b)*mesh.magSf()));
 
     volScalarField mgb("mgb", mag(n));
-    volScalarField ddtb("ddtb",fvc::ddt(b));
 
     dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL);
 
     {
         volScalarField bc(b*c);
 
-        dMgb += 1.0e-8*
+        dMgb += 1.0e-3*
             (bc*mgb)().weightedAverage(mesh.V())
            /(bc.weightedAverage(mesh.V()) + SMALL);
     }
@@ -63,7 +54,6 @@ if (ign.ignited())
     surfaceScalarField nf("nf", mesh.Sf() & nfVec);
     n /= mgb;
 
-
     #include "StCorr.H"
 
     // Calculate turbulent flame speed flux
@@ -74,6 +64,8 @@ if (ign.ignited())
 
     Db = flameWrinkling->Db();
 
+    // Create b equation
+    // ~~~~~~~~~~~~~~~~~
     fvScalarMatrix bEqn
     (
         betav*fvm::ddt(rho, b)
@@ -83,7 +75,6 @@ if (ign.ignited())
       - fvm::laplacian(Db, b)
      ==
         betav*fvOptions(rho, b)
-        // + betav*rhou*St*StCorr*mag(fvc::grad(b))
     );
 
 
@@ -91,14 +82,17 @@ if (ign.ignited())
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     #include "ignite.H"
 
+    // Solve for b
+    // ~~~~~~~~~~~
     bEqn.relax();
+
     fvOptions.constrain(bEqn);
+
     bEqn.solve();
 
     fvOptions.correct(b);
 
-    Info<< "min(b) = " << min(b).value() << nl
-        << "max(b) = " << max(b).value() << endl;
+    Info<< "min(b) = " << min(b).value() << endl;
 
     if (composition.contains("ft"))
     {
@@ -115,14 +109,8 @@ if (ign.ignited())
             << endl;
     }
 
-
-    // Correct the flame-wrinkling, smoothing
+    // Correct the flame-wrinkling
     flameWrinkling->correct(mvConvection);
-    St = Xp*Xi*Su*combustFlag;
 
-    if (PDRProperties.getOrDefault("smoothSt", false))
-    {
-        const scalar coeff = PDRProperties.get<scalar>("smoothStCoeff");
-        fvc::smooth(St, coeff);
-    }
+    St = Xi*Su;
 }

applications/solvers/combustion/PDRFoam/createFieldRefs.H

@@ -1,2 +1 @@
 const volScalarField& psi = thermo.psi();
-scalar pMin = thermo.getOrDefault<scalar>("pMin", 0);

applications/solvers/combustion/PDRFoam/createFields.H

@@ -1,303 +1,218 @@
-    Info<< "Reading thermophysical properties\n" << endl;
+Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<psiuReactionThermo> pThermo
+autoPtr<psiuReactionThermo> pThermo
+(
+    psiuReactionThermo::New(mesh)
+);
+psiuReactionThermo& thermo = pThermo();
+thermo.validate(args.executable(), "ha", "ea");
+
+basicSpecieMixture& composition = thermo.composition();
+
+volScalarField rho
+(
+    IOobject
     (
-        psiuReactionThermo::New(mesh)
-    );
-    psiuReactionThermo& thermo = pThermo();
-    thermo.validate(args.executable(), "ha", "ea");
-    basicMultiComponentMixture& composition = thermo.composition();
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
 
-    volScalarField rho
+volScalarField& p = thermo.p();
+
+volScalarField& b = composition.Y("b");
+Info<< "min(b) = " << min(b).value() << endl;
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
     (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo.rho()
-    );
-    volScalarField rhou
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+mesh.setFluxRequired(p.name());
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::RASModel> turbulence
+(
+    compressible::New<compressible::RASModel>
     (
-        IOobject
-        (
-            "rhou",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo.rhou()
-    );
-    volScalarField muu
-    (
-        IOobject
-        (
-            "muu",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo.muu()
-    );
-
-    volScalarField& p = thermo.p();
-
-    volScalarField& b = composition.Y("b");
-
-    volScalarField combustFlag
-    (
-        IOobject
-        (
-            "combustFlag",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading fields betav, betai\n" << endl;
-    volScalarField betav
-    (
-        IOobject
-        (
-            "betav",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    volSymmTensorField betai
-    (
-        IOobject
-        (
-            "betai",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field Bv\n" << endl;
-    volSymmTensorField Bv
-    (
-        IOobject
-        (
-            "Bv",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading fields Nv,nsv\n" << endl;
-    volScalarField Nv
-    (
-        IOobject
-        (
-            "Nv",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    volSymmTensorField nsv
-    (
-        IOobject
-        (
-            "nsv",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading fields Blong\n" << endl;
-    volSymmTensorField Blong
-    (
-        IOobject
-        (
-            "Blong",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "\nCreating field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "compressibleCreatePhi.H"
-    mesh.setFluxRequired(p.name());
-
-    Info<< "Creating turbulence model\n" << endl;
-    autoPtr<compressible::RASModel> turbulence
-    (
-        compressible::New<compressible::RASModel>
-        (
-            rho,
-            U,
-            phi,
-            thermo
-        )
-    );
-
-    #include "createDpdt.H"
-
-    #include "createK.H"
-
-    Info<< "Creating the unstrained laminar flame speed\n" << endl;
-    autoPtr<laminarFlameSpeed> unstrainedLaminarFlameSpeed
-    (
-        laminarFlameSpeed::New(thermo)
-    );
-
-    Info<< "Reading strained laminar flame speed field Su\n" << endl;
-    volScalarField Su
-    (
-        IOobject
-        (
-            "Su",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Reading field Lobs\n" << endl;
-    volScalarField Lobs
-    (
-        IOobject
-        (
-            "Lobs",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field CT\n" << endl;
-    volSymmTensorField CT
-    (
-        IOobject
-        (
-            "CT",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    IOdictionary PDRProperties
-    (
-        IOobject
-        (
-            "PDRProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-    //- Create the drag model
-    autoPtr<PDRDragModel> drag = PDRDragModel::New
-    (
-        PDRProperties,
-        turbulence,
         rho,
         U,
-        phi
-    );
+        phi,
+        thermo
+    )
+);
 
-    //- Create the flame-wrinkling model
-    autoPtr<XiModel> flameWrinkling = XiModel::New
+
+Info<< "Creating the unstrained laminar flame speed\n" << endl;
+autoPtr<laminarFlameSpeed> unstrainedLaminarFlameSpeed
+(
+    laminarFlameSpeed::New(thermo)
+);
+
+
+Info<< "Reading strained laminar flame speed field Su\n" << endl;
+volScalarField Su
+(
+    IOobject
     (
-        PDRProperties,
-        thermo,
-        turbulence,
-        Su,
-        rho,
-        b,
-        phi
-    );
+        "Su",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
 
-    Info<< "Calculating turbulent flame speed field St\n" << endl;
-    volScalarField St
+Info<< "Reading field betav\n" << endl;
+volScalarField betav
+(
+    IOobject
     (
-        IOobject
-        (
-            "St",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        flameWrinkling->Xi()*Su
-    );
+        "betav",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
 
-    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-    if (composition.contains("ft"))
-    {
-        fields.add(composition.Y("ft"));
-    }
-
-    fields.add(b);
-    fields.add(thermo.he());
-    fields.add(thermo.heu());
-    flameWrinkling->addXi(fields);
-
-    #include "createMRF.H"
-    #include "createFvOptions.H"
-
-    volScalarField c
+Info<< "Reading field Lobs\n" << endl;
+volScalarField Lobs
+(
+    IOobject
     (
-        IOobject
-        (
-            "c",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        scalar(1) - b
-    );
+        "Lobs",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field CT\n" << endl;
+volSymmTensorField CT
+(
+    IOobject
+    (
+        "CT",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field Nv\n" << endl;
+volScalarField Nv
+(
+    IOobject
+    (
+        "Nv",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field nsv\n" << endl;
+volSymmTensorField nsv
+(
+    IOobject
+    (
+        "nsv",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+IOdictionary PDRProperties
+(
+    IOobject
+    (
+        "PDRProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+//- Create the drag model
+autoPtr<PDRDragModel> drag = PDRDragModel::New
+(
+    PDRProperties,
+    *turbulence,
+    rho,
+    U,
+    phi
+);
+
+//- Create the flame-wrinkling model
+autoPtr<XiModel> flameWrinkling = XiModel::New
+(
+    PDRProperties,
+    thermo,
+    *turbulence,
+    Su,
+    rho,
+    b,
+    phi
+);
+
+Info<< "Calculating turbulent flame speed field St\n" << endl;
+volScalarField St
+(
+    IOobject
+    (
+        "St",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    flameWrinkling->Xi()*Su
+);
+
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+if (composition.contains("ft"))
+{
+    fields.add(composition.Y("ft"));
+}
+
+fields.add(b);
+fields.add(thermo.he());
+fields.add(thermo.heu());
+flameWrinkling->addXi(fields);
+
+#include "createDpdt.H"
+
+#include "createK.H"
+
+#include "createMRF.H"
+#include "createFvOptions.H"

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -75,28 +75,32 @@ Foam::laminarFlameSpeedModels::SCOPE::SCOPE
     (
         dictionary
         (
-            IFstream(dict.get<fileName>("fuelFile"))()
-        ).subDict(typeName + "Coeffs")
+            IFstream
+            (
+                dict.get<fileName>("fuelFile")
+            )()
+        ).optionalSubDict(typeName + "Coeffs")
+    ),
+    LFL_
+    (
+        coeffsDict_.getCompat<scalar>
+        (
+            "lowerFlammabilityLimit",
+            {{"lowerFlamabilityLimit", 1712}}
+        )
+    ),
+    UFL_
+    (
+        coeffsDict_.getCompat<scalar>
+        (
+            "upperFlammabilityLimit",
+            {{"upperFlamabilityLimit", 1712}}
+        )
     ),
-    LFL_(coeffsDict_.get<scalar>("lowerFlamabilityLimit")),
-    UFL_(coeffsDict_.get<scalar>("upperFlamabilityLimit")),
     SuPolyL_(coeffsDict_.subDict("lowerSuPolynomial")),
     SuPolyU_(coeffsDict_.subDict("upperSuPolynomial")),
-    Texp0_(coeffsDict_.get<scalar>("Texp0")),
-    Texp1_(coeffsDict_.get<scalar>("Texp1")),
-    Texp2_(coeffsDict_.get<scalar>("Texp2")),
-    Texp3_(coeffsDict_.get<scalar>("Texp3")),
-    Texp4_(coeffsDict_.get<scalar>("Texp4")),
-    Texp5_(coeffsDict_.get<scalar>("Texp5")),
-    Texp6_(coeffsDict_.get<scalar>("Texp6")),
-    pexp0_(coeffsDict_.get<scalar>("pexp0")),
-    pexp1_(coeffsDict_.get<scalar>("pexp1")),
-    pexp2_(coeffsDict_.get<scalar>("pexp2")),
-    pexp3_(coeffsDict_.get<scalar>("pexp3")),
-    pexp4_(coeffsDict_.get<scalar>("pexp4")),
-    pexp5_(coeffsDict_.get<scalar>("pexp5")),
-    pexp6_(coeffsDict_.get<scalar>("pexp6")),
-    CIn_(coeffsDict_.getOrDefault<scalar>("CIn", 0)),
+    Texp_(coeffsDict_.get<scalar>("Texp")),
+    pexp_(coeffsDict_.get<scalar>("pexp")),
     MaPolyL_(coeffsDict_.subDict("lowerMaPolynomial")),
     MaPolyU_(coeffsDict_.subDict("upperMaPolynomial"))
 {
@@ -154,19 +158,19 @@ inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::SuRef
 {
     if (phi < LFL_ || phi > UFL_)
     {
-        // Return 0 beyond the flamibility limits
+        // Return 0 beyond the flammability limits
         return scalar(0);
     }
     else if (phi < SuPolyL_.ll)
     {
         // Use linear interpolation between the low end of the
-        // lower polynomial and the lower flamibility limit
+        // lower polynomial and the lower flammability limit
         return SuPolyL_.llv*(phi - LFL_)/(SuPolyL_.ll - LFL_);
     }
     else if (phi > SuPolyU_.ul)
     {
         // Use linear interpolation between the upper end of the
-        // upper polynomial and the upper flamibility limit
+        // upper polynomial and the upper flammability limit
         return SuPolyU_.ulv*(UFL_ - phi)/(UFL_ - SuPolyU_.ul);
     }
     else if (phi < SuPolyL_.lu)
@@ -181,7 +185,7 @@ inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::SuRef
     }
     else
     {
-        FatalErrorIn("laminarFlameSpeedModels::SCOPE::SuRef(scalar phi)")
+        FatalErrorInFunction
             << "phi = " << phi
             << " cannot be handled by SCOPE function with the "
                "given coefficients"
@@ -191,41 +195,35 @@ inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::SuRef
     }
 }
 
+
 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::Ma
 (
     scalar phi
 ) const
 {
-    if (phi < LFL_ || phi > UFL_)
+    if (phi < MaPolyL_.ll)
     {
-        // Return 0 beyond the flamibility limits
-        return scalar(0);
-    }
-    else if (phi < MaPolyL_.ll)
-    {
-        // Use linear interpolation between the low end of the
-        // lower polynomial and the lower flamibility limit
-        return MaPolyL_.llv*(phi - LFL_)/(MaPolyL_.ll - LFL_);
+        // Beyond the lower limit assume Ma is constant
+        return MaPolyL_.llv;
     }
     else if (phi > MaPolyU_.ul)
     {
-        // Use linear interpolation between the upper end of the
-        // upper polynomial and the upper flamibility limit
-        return MaPolyU_.ulv*(UFL_ - phi)/(UFL_ - MaPolyU_.ul);
+        // Beyond the upper limit assume Ma is constant
+        return MaPolyU_.ulv;
     }
-    else if (phi < MaPolyL_.lu)
+    else if (phi < SuPolyL_.lu)
     {
         // Evaluate the lower polynomial
         return polyPhi(phi, MaPolyL_);
     }
-    else if (phi > MaPolyU_.lu)
+    else if (phi > SuPolyU_.lu)
     {
         // Evaluate the upper polynomial
         return polyPhi(phi, MaPolyU_);
     }
     else
     {
-        FatalErrorIn("laminarFlameSpeedModels::SCOPE::Ma(scalar phi)")
+        FatalErrorInFunction
             << "phi = " << phi
             << " cannot be handled by SCOPE function with the "
                "given coefficients"
@@ -245,31 +243,8 @@ inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
 {
     static const scalar Tref = 300.0;
     static const scalar pRef = 1.013e5;
-    scalar phi1=phi-1.0;
 
-    scalar Texp_=
-        Texp0_+phi1*
-        (
-            Texp1_
-            +Texp2_*phi1+
-            Texp3_*pow(phi1, 2) +
-            Texp4_*pow(phi1, 3) +
-            Texp5_*pow(phi1, 4) +
-            Texp6_*pow(phi1, 5)
-        );
-
-    scalar pexp =
-        pexp0_ + phi1
-       *(
-            pexp1_
-          + pexp2_*phi1
-          + pexp3_*pow(phi1, 2)
-          + pexp4_*pow(phi1, 3)
-          + pexp5_*pow(phi1, 4)
-          + pexp6_*pow(phi1, 5)
-        );
-
-    return SuRef(phi)*pow((Tu/Tref), Texp_)*pow((p/pRef), pexp);
+    return SuRef(phi)*pow((Tu/Tref), Texp_)*pow((p/pRef), pexp_);
 }
 
 
@@ -304,9 +279,11 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi);
     }
 
-    forAll(Su0.boundaryField(), patchi)
+    volScalarField::Boundary& Su0Bf = Su0.boundaryFieldRef();
+
+    forAll(Su0Bf, patchi)
     {
-        scalarField& Su0p = Su0.boundaryFieldRef()[patchi];
+        scalarField& Su0p = Su0Bf[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
         const scalarField& Tup = Tu.boundaryField()[patchi];
 
@@ -351,9 +328,11 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
         Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi[celli]);
     }
 
-    forAll(Su0.boundaryField(), patchi)
+    volScalarField::Boundary& Su0Bf = Su0.boundaryFieldRef();
+
+    forAll(Su0Bf, patchi)
     {
-        scalarField& Su0p = Su0.boundaryFieldRef()[patchi];
+        scalarField& Su0p = Su0Bf[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
         const scalarField& Tup = Tu.boundaryField()[patchi];
         const scalarField& phip = phi.boundaryField()[patchi];
@@ -373,6 +352,7 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
     return tSu0;
 }
 
+
 Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
 (
     const volScalarField& phi
@@ -402,11 +382,11 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
         ma[celli] = Ma(phi[celli]);
     }
 
+    volScalarField::Boundary& maBf = ma.boundaryFieldRef();
 
-
-    forAll(ma.boundaryField(), patchi)
+    forAll(maBf, patchi)
     {
-        scalarField& map = ma.boundaryFieldRef()[patchi];
+        scalarField& map = maBf[patchi];
         const scalarField& phip = phi.boundaryField()[patchi];
 
         forAll(map, facei)
@@ -418,11 +398,6 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
     return tMa;
 }
 
-Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::CIn() const
-{
-    return  CIn_ ;
-}
-
 
 Foam::tmp<Foam::volScalarField>
 Foam::laminarFlameSpeedModels::SCOPE::Ma() const
@@ -435,8 +410,7 @@ Foam::laminarFlameSpeedModels::SCOPE::Ma() const
         (
             dimensionedScalar
             (
-                "stoichiometricAirFuelMassRatio",
-                psiuReactionThermo_
+                "stoichiometricAirFuelMassRatio", dimless, psiuReactionThermo_
             )*ft/(scalar(1) - ft)
         );
     }
@@ -477,8 +451,7 @@ Foam::laminarFlameSpeedModels::SCOPE::operator()() const
             psiuReactionThermo_.Tu(),
             dimensionedScalar
             (
-                "stoichiometricAirFuelMassRatio",
-                psiuReactionThermo_
+                "stoichiometricAirFuelMassRatio", dimless, psiuReactionThermo_
             )*ft/(scalar(1) - ft)
         );
     }

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H

@@ -123,50 +123,11 @@ class SCOPE
         //- Upper Su polynomial
         polynomial SuPolyU_;
 
-        //- Temperature correction exponent0
-        scalar Texp0_;
+        //- Temperature correction exponent
+        scalar Texp_;
 
-        //- Temperature correction exponent1
-        scalar Texp1_;
-
-        //- Temperature correction exponent2
-        scalar Texp2_;
-
-        //- Temperature correction exponent3
-        scalar Texp3_;
-
-        //- Temperature correction exponent4
-        scalar Texp4_;
-
-        //- Temperature correction exponent5
-        scalar Texp5_;
-
-        //- Temperature correction exponent6
-        scalar Texp6_;
-
-        //- Pressure correction exponent0
-        scalar pexp0_;
-
-        //- Pressure correction exponent1
-        scalar pexp1_;
-
-        //- pressure correction exponent2
-        scalar pexp2_;
-
-        //- pressure correction exponent3
-        scalar pexp3_;
-
-        //- Pressure correction exponent4
-        scalar pexp4_;
-
-        //- Pressure correction exponent5
-        scalar pexp5_;
-
-        //- pressure correction exponent6
-        scalar pexp6_;
-
-        //- Additional coefficient
-        scalar CIn_;
+        //- Pressure correction exponent
+        scalar pexp_;
 
         //- Lower Ma polynomial
         polynomial MaPolyL_;
@@ -175,25 +136,25 @@ class SCOPE
         polynomial MaPolyU_;
 
 
-    // Private Member Functions
+    // Private member functions
 
         //- Polynomial evaluated from the given equivalence ratio
-        //- and polynomial coefficients
+        //  and polynomial coefficients
         static inline scalar polyPhi(scalar phi, const polynomial& a);
 
         //- Laminar flame speed evaluated from the given equivalence ratio
-        //- at the reference temperature and pressure
+        //  at the reference temperature and pressure
         inline scalar SuRef(scalar phi) const;
 
         //- Markstein evaluated from the given equivalence ratio
         inline scalar Ma(scalar phi) const;
 
         //- Laminar flame speed evaluated from the given equivalence ratio
-        //- corrected for temperature and pressure dependence
+        //  corrected for temperature and pressure dependence
         inline scalar Su0pTphi(scalar p, scalar Tu, scalar phi) const;
 
         //- Laminar flame speed evaluated from the given uniform
-        //- equivalence ratio corrected for temperature and pressure dependence
+        //  equivalence ratio corrected for temperature and pressure dependence
         tmp<volScalarField> Su0pTphi
         (
             const volScalarField& p,
@@ -202,7 +163,7 @@ class SCOPE
         ) const;
 
         //- Laminar flame speed evaluated from the given equivalence ratio
-        //- distribution corrected for temperature and pressure dependence
+        //  distribution corrected for temperature and pressure dependence
         tmp<volScalarField> Su0pTphi
         (
             const volScalarField& p,
@@ -211,14 +172,13 @@ class SCOPE
         ) const;
 
         //- Return the Markstein number
-        //- evaluated from the given equivalence ratio
+        //  evaluated from the given equivalence ratio
         tmp<volScalarField> Ma(const volScalarField& phi) const;
 
-        //- No copy construct
-        SCOPE(const SCOPE&) = delete;
+        //- Construct as copy (not implemented)
+        SCOPE(const SCOPE&);
 
-        //- No copy assign
-        void operator=(const SCOPE&) = delete;
+        void operator=(const SCOPE&);
 
 
 public:
@@ -231,8 +191,8 @@ public:
         //- Construct from dictionary and psiuReactionThermo
         SCOPE
         (
-            const dictionary& dict,
-            const psiuReactionThermo& ct
+            const dictionary&,
+            const psiuReactionThermo&
         );
 
 
@@ -245,9 +205,6 @@ public:
         //- Return the Markstein number
         tmp<volScalarField> Ma() const;
 
-        //- Return CIn
-        scalar CIn() const;
-
         //- Return the laminar flame speed [m/s]
         tmp<volScalarField> operator()() const;
 };

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -1,7 +1,6 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-
 volVectorField HbyA(constrainHbyA(invA & UEqn.H(), U, p));
 
 if (pimple.transonic())
@@ -9,7 +8,12 @@ if (pimple.transonic())
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)*(fvc::flux(HbyA))
+        fvc::interpolate(psi)
+       *(
+            fvc::flux(HbyA)
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+           /fvc::interpolate(rho)
+        )
     );
 
     while (pimple.correctNonOrthogonal())
@@ -36,7 +40,10 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::flux(rho*HbyA)
+        (
+            fvc::flux(rho*HbyA)
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+        )
     );
 
     while (pimple.correctNonOrthogonal())
@@ -59,8 +66,6 @@ else
     }
 }
 
-p.max(pMin);
-
 #include "rhoEqn.H"
 #include "continuityErrs.H"
 

applications/solvers/combustion/PDRFoam/readTimeControls.H

@@ -1,47 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Global
-    readTimeControls
-
-Description
-    Read the control parameters used by setDeltaT
-
-\*---------------------------------------------------------------------------*/
-
-const bool adjustTimeStep =
-    runTime.controlDict().getOrDefault<bool>("adjustTimeStep", false);
-
-scalar maxCo = runTime.controlDict().getOrDefault<scalar>("maxCo", 1.0);
-
-scalar maxDeltaT =
-    runTime.controlDict().getOrDefault<scalar>("maxDeltaT", GREAT);
-
-scalar maxDeltaTRatio =
-    runTime.controlDict().getOrDefault<scalar>("maxDeltaTRatio", 1.2);
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/rhoEqn.H

@@ -1,3 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
 {
     solve(betav*fvm::ddt(rho) + fvc::div(phi));
 }
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/setDeltaT.H

@@ -28,7 +28,7 @@ Global
 
 Description
     Reset the timestep to maintain a constant maximum courant Number.
-    Reduction of time-step is immediate but increase is damped to avoid
+    Reduction of time-step is imediate but increase is damped to avoid
     unstable oscillations.
 
 \*---------------------------------------------------------------------------*/
@@ -36,8 +36,7 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact = maxCo/(CoNum + StCoNum + SMALL);
-    scalar deltaTFact =
-        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), maxDeltaTRatio);
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (

applications/solvers/combustion/XiFoam/XiDyMFoam/XiDyMFoam.C

@@ -36,7 +36,7 @@ Description
 
     Combusting RANS code using the b-Xi two-equation model.
     Xi may be obtained by either the solution of the Xi transport
-    equation or from an algebraic expression.  Both approaches are
+    equation or from an algebraic exression.  Both approaches are
     based on Gulder's flame speed correlation which has been shown
     to be appropriate by comparison with the results from the
     spectral model.

applications/solvers/combustion/XiFoam/XiDyMFoam/createControls.H

@@ -2,10 +2,10 @@
 
 bool correctPhi
 (
-    pimple.dict().getOrDefault("correctPhi", true)
+    pimple.dict().lookupOrDefault("correctPhi", true)
 );
 
 bool checkMeshCourantNo
 (
-    pimple.dict().getOrDefault("checkMeshCourantNo", false)
+    pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
 );

applications/solvers/combustion/XiFoam/XiDyMFoam/readControls.H

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

applications/solvers/combustion/XiFoam/XiEngineFoam/XiEngineFoam.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +31,7 @@ Description
 
     Combusting RANS code using the b-Xi two-equation model.
     Xi may be obtained by either the solution of the Xi transport
-    equation or from an algebraic expression.  Both approaches are
+    equation or from an algebraic exression.  Both approaches are
     based on Gulder's flame speed correlation which has been shown
     to be appropriate by comparison with the results from the
     spectral model.
@@ -145,7 +144,9 @@ int main(int argc, char *argv[])
 
         runTime.write();
 
-        runTime.printExecutionTime(Info);
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
     }
 
     Info<< "End\n" << endl;

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -35,7 +35,7 @@ Description
 
     Combusting RANS code using the b-Xi two-equation model.
     Xi may be obtained by either the solution of the Xi transport
-    equation or from an algebraic expression.  Both approaches are
+    equation or from an algebraic exression.  Both approaches are
     based on Gulder's flame speed correlation which has been shown
     to be appropriate by comparison with the results from the
     spectral model.

applications/solvers/combustion/chemFoam/Make/options

@@ -6,13 +6,11 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lfluidThermophysicalModels \
     -lreactionThermophysicalModels \
-    -lchemistryModel \
-    -lthermophysicalProperties
+    -lchemistryModel

applications/solvers/combustion/chemFoam/chemFoam.C

@@ -93,8 +93,8 @@ int main(int argc, char *argv[])
         runTime.printExecutionTime(Info);
     }
 
-    Info<< "Number of steps = " << runTime.timeIndex() << nl;
-    Info<< "End\n" << endl;
+    Info << "Number of steps = " << runTime.timeIndex() << endl;
+    Info << "End" << nl << endl;
 
     return 0;
 }

applications/solvers/combustion/chemFoam/readInitialConditions.H

@@ -37,7 +37,11 @@
     {
         forAll(Y, i)
         {
-            fractions.readIfPresent(Y[i].name(), X0[i]);
+            const word& name = Y[i].name();
+            if (fractions.found(name))
+            {
+                X0[i] = fractions.get<scalar>(name);
+            }
         }
 
         scalar mw = 0.0;
@@ -57,7 +61,11 @@
     {
         forAll(Y, i)
         {
-            fractions.readIfPresent(Y[i].name(), Y0[i]);
+            const word& name = Y[i].name();
+            if (fractions.found(name))
+            {
+                Y0[i] = fractions.get<scalar>(name);
+            }
         }
 
         scalar invW = 0.0;