diff --git a/.clang-format b/.clang-format new file mode 100644 index 00000000..a37fea63 --- /dev/null +++ b/.clang-format @@ -0,0 +1,66 @@ +--- +Language: Cpp +BasedOnStyle: LLVM + +ColumnLimit: 80 +IndentWidth: 4 +ContinuationIndentWidth: 4 +TabWidth: 4 +UseTab: Never + +# OpenFOAM-style brace placement. +BreakBeforeBraces: Custom +BraceWrapping: + AfterClass: true + AfterControlStatement: Always + AfterEnum: true + AfterFunction: true + AfterNamespace: true + AfterStruct: true + AfterUnion: true + BeforeCatch: true + BeforeElse: true + SplitEmptyFunction: false + SplitEmptyRecord: false + SplitEmptyNamespace: false + +AllowShortBlocksOnASingleLine: Never +AllowShortCaseLabelsOnASingleLine: false +AllowShortEnumsOnASingleLine: false +AllowShortFunctionsOnASingleLine: None +AllowShortIfStatementsOnASingleLine: Never +AllowShortLoopsOnASingleLine: false + +# if (...), for (...), while (...), switch (...) +# but: +# forAll(...), not forAll (...) +SpaceBeforeParens: Custom +SpaceBeforeParensOptions: + AfterControlStatements: true + AfterForeachMacros: false + AfterFunctionDeclarationName: false + AfterFunctionDefinitionName: false + AfterIfMacros: false + AfterOverloadedOperator: false + AfterPlacementOperator: false + +ForEachMacros: + - forAll + - forAllReverse + - forAllIter + - forAllConstIter + - forAllIters + - forAllConstIters + +DerivePointerAlignment: false +PointerAlignment: Left +ReferenceAlignment: Left + +SortIncludes: + Enabled: false +ReflowComments: Never +FixNamespaceComments: false + +BreakBeforeBinaryOperators: NonAssignment +AlignOperands: AlignAfterOperator +... diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index a8b0756e..6bc2e9c0 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -21,21 +21,21 @@ jobs: runs-on: ubuntu-latest container: - image: docker.io/openfoam/openfoam10-paraview510 + image: microfluidica/openfoam:13 options: --user root steps: - name: Checkout AdditiveFOAM uses: actions/checkout@v2 - name: Build AdditiveFOAM run: | - . /opt/openfoam10/etc/bashrc + . /opt/openfoam13/etc/bashrc || true + test -n "$WM_PROJECT_DIR" ./Allwmake - name: Test AdditiveFOAM run: | - . /opt/openfoam10/etc/bashrc + . /opt/openfoam13/etc/bashrc || true cp -r tutorials/AMB2018-02-B userCase cd userCase - # FIXME: use built-in "additiveFoam" smaller case when created blockMesh decomposePar mpirun -n 6 --oversubscribe --allow-run-as-root additiveFoam -parallel diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..18d24c1f --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,67 @@ +minimum_pre_commit_version: "3.7.0" + +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v6.0.0 + hooks: + - id: trailing-whitespace + files: &openfoam_cpp_files '\.(C|H|I|cc|cpp|cxx|hh|hpp|hxx)$' + exclude: &openfoam_excludes '(^|/)(lnInclude|ThirdParty|platforms|build|processor[0-9]+|postProcessing)/' + + - id: end-of-file-fixer + files: *openfoam_cpp_files + exclude: *openfoam_excludes + + - id: mixed-line-ending + args: [--fix=lf] + files: *openfoam_cpp_files + exclude: *openfoam_excludes + + - id: check-merge-conflict + + - id: check-yaml + + - repo: local + hooks: + - id: openfoam-no-tabs + name: "OpenFOAM style: no tabs" + language: pygrep + entry: '\t' + files: *openfoam_cpp_files + exclude: *openfoam_excludes + + # Manual for now because fixing long C++ lines usually requires judgment. + # Remove "stages: [manual]" once the existing codebase is clean. + - id: openfoam-max-80-columns + name: "OpenFOAM style: max 80 characters" + language: pygrep + entry: '^.{81,}$' + files: *openfoam_cpp_files + exclude: *openfoam_excludes + stages: [manual] + + - id: openfoam-space-after-control-keyword + name: "OpenFOAM style: use if (...), for (...), while (...), switch (...)" + language: pygrep + entry: '^(?!\s*#).*\b(if|for|while|switch)\(' + files: *openfoam_cpp_files + exclude: *openfoam_excludes + + - id: openfoam-forall-no-space + name: "OpenFOAM style: use forAll(...), not forAll (...)" + language: pygrep + entry: '\b(forAll|forAllReverse|forAllIter|forAllConstIter|forAllIters|forAllConstIters)\s+\(' + files: *openfoam_cpp_files + exclude: *openfoam_excludes + + # Manual only: clang-format is useful as a helper, but it cannot exactly + # represent OpenFOAM stream/math/operator spacing conventions. + - repo: https://github.com/pre-commit/mirrors-clang-format + rev: v22.1.5 + hooks: + - id: clang-format + name: "clang-format: approximate OpenFOAM style, requires manual changes" + args: [--style=file] + files: *openfoam_cpp_files + exclude: *openfoam_excludes + stages: [manual] diff --git a/README.md b/README.md index 85ae7fc7..ae69d800 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ The documentation for `AdditiveFOAM` is hosted on [GitHub Pages](https://ornl.gi | [tutorials](tutorials) | Tutorial cases | ## Installation and Dependencies -[![OpenFOAM-10](https://img.shields.io/badge/OpenFOAM-10-blue.svg)](https://github.com/OpenFOAM/OpenFOAM-10) +[![OpenFOAM-13](https://img.shields.io/badge/OpenFOAM-13-blue.svg)](https://github.com/OpenFOAM/OpenFOAM-13) AdditiveFOAM is built on source code released by the OpenFOAM Foundation [openfoam.org](https://openfoam.org/), which is available in public [OpenFOAM repositories](https://github.com/OpenFOAM). diff --git a/applications/solvers/additiveFoam/Allwmake b/applications/solvers/additiveFoam/Allwmake index e54e7e16..9735535b 100755 --- a/applications/solvers/additiveFoam/Allwmake +++ b/applications/solvers/additiveFoam/Allwmake @@ -28,6 +28,7 @@ export ADDITIVEFOAM_BUILD_FLAGS="-DGIT_MODULE_ENABLED=1" #------------------------------------------------------------------------------ # Build libraries and solver wmake $targetType functionObjects +wmake $targetType utilities wmake $targetType movingHeatSource wmake $targetType diff --git a/applications/solvers/additiveFoam/Make/options b/applications/solvers/additiveFoam/Make/options index c61a6445..2e81c929 100644 --- a/applications/solvers/additiveFoam/Make/options +++ b/applications/solvers/additiveFoam/Make/options @@ -2,13 +2,16 @@ EXE_INC = \ $(ADDITIVEFOAM_BUILD_FLAGS) \ -I. \ -ImovingHeatSource/lnInclude \ + -Iutilities/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/randomProcesses/lnInclude EXE_LIBS = \ -L$(FOAM_USER_LIBBIN) \ -lmovingBeamModels \ + -ladditiveFoamUtilities \ -lfiniteVolume \ -lmeshTools \ -lrandomProcesses diff --git a/applications/solvers/additiveFoam/additiveFoam.C b/applications/solvers/additiveFoam/additiveFoam.C index 20d70b1a..a12dfadc 100644 --- a/applications/solvers/additiveFoam/additiveFoam.C +++ b/applications/solvers/additiveFoam/additiveFoam.C @@ -5,7 +5,7 @@ \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2023 Oak Ridge National Laboratory + Copyright (C) 2023 Oak Ridge National Laboratory ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -29,34 +29,57 @@ Application Description A transient heat transfer and fluid flow solver for additive manufacturing simulations. - + \*---------------------------------------------------------------------------*/ -#include "additiveFoamInfo.H" -#include "fvCFD.H" + +#include "argList.H" +#include "timeSelector.H" +#include "zeroGradientFvPatchFields.H" +#include "IFstream.H" +#include "uniformDimensionedFields.H" +#include "pressureReference.H" +#include "findRefCell.H" + +#include "fvmDiv.H" +#include "fvmDdt.H" +#include "fvcSurfaceIntegrate.H" +#include "fvcVolumeIntegrate.H" +#include "fvmLaplacian.H" +#include "constrainPressure.H" +#include "constrainHbyA.H" +#include "adjustPhi.H" #include "pimpleControl.H" -#include "graph.H" +#include "fvCorrectPhi.H" #include "Polynomial.H" -#include "interpolateXY/interpolateXY.H" -#include "movingHeatSourceModel.H" + #include "EulerDdtScheme.H" #include "CrankNicolsonDdtScheme.H" +#include "additiveFoamInfo.H" +#include "movingHeatSourceModel.H" +#include "graph.H" +#include "interpolateXY.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) -{ - #include "postProcess.H" +{ + using namespace Foam; + + #include "postProcess.H" #include "setRootCase.H" + AdditiveFoamInfo::write(); + #include "createTime.H" #include "createMesh.H" - #include "createControl.H" + #include "createDyMControls.H" #include "createFields.H" - #include "createTimeControls.H" #include "initContinuityErrs.H" - + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + // Initialize time-stepping controls scalar DiNum = 0.0; scalar alphaCoNum = 0.0; movingHeatSourceModel sources(mesh); @@ -67,26 +90,33 @@ int main(int argc, char *argv[]) { #include "updateProperties.H" - #include "readTimeControls.H" + #include "readDyMControls.H" #include "CourantNo.H" #include "setDeltaT.H" sources.update(); - + + mesh.update(); + runTime++; - Info<< "Time = " << runTime.timeName() << nl << endl; + Info<< "Time = " << runTime.name() << nl << endl; #include "solutionControls.H" - - while (pimple.loop() && fluidInDomain) + + while (pimple.loop()) { - #include "pU/UEqn.H" - #include "pU/pEqn.H" + #include "moveMesh.H" + + if (fluidInDomain) + { + #include "pU/UEqn.H" + #include "pU/pEqn.H" + } } #include "thermo/TEqn.H" - + runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" diff --git a/applications/solvers/additiveFoam/additiveFoamInfo.H b/applications/solvers/additiveFoam/additiveFoamInfo.H index 88cea642..8ea78586 100644 --- a/applications/solvers/additiveFoam/additiveFoamInfo.H +++ b/applications/solvers/additiveFoam/additiveFoamInfo.H @@ -1,6 +1,6 @@ /*---------------------------------------------------------------------------*\ ------------------------------------------------------------------------------- - Copyright (C) 2023 Oak Ridge National Laboratory + Copyright (C) 2023 Oak Ridge National Laboratory ------------------------------------------------------------------------------- Class @@ -21,7 +21,7 @@ Description #endif // Static version -#define ADDITIVEFOAM_VERSION "1.2.0-dev" +#define ADDITIVEFOAM_VERSION "2.0.0-dev" namespace Foam { @@ -60,4 +60,3 @@ public: #endif // ************************************************************************* // - diff --git a/applications/solvers/additiveFoam/createFields.H b/applications/solvers/additiveFoam/createFields.H index a197af73..4a45b0ba 100644 --- a/applications/solvers/additiveFoam/createFields.H +++ b/applications/solvers/additiveFoam/createFields.H @@ -4,7 +4,7 @@ volScalarField T IOobject ( "T", - runTime.timeName(), + runTime.name(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE @@ -18,7 +18,7 @@ volScalarField alpha3 IOobject ( IOobject::groupName("alpha", "powder"), - runTime.timeName(), + runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE @@ -33,7 +33,7 @@ volScalarField alpha1 IOobject ( IOobject::groupName("alpha", "solid"), - runTime.timeName(), + runTime.name(), mesh, IOobject::NO_READ, IOobject::NO_WRITE @@ -45,7 +45,8 @@ volScalarField alpha1 IFstream thermoFile ( - runTime.rootPath()/runTime.globalCaseName()/runTime.constant()/"thermoPath" + runTime.rootPath()/runTime.globalCaseName() + /runTime.constant()/"thermoPath" ); graph thermo @@ -71,10 +72,10 @@ dimensionedScalar Tsol ); // set solid fraction field consistent with temperature -forAll(mesh.cells(), cellI) +forAll(mesh.cells(), celli) { - scalar alpha1_ = interpolateXY(T[cellI], thermo.x(), thermo.y()); - alpha1[cellI] = min(max(alpha1_, 0.0), 1.0); + scalar alpha1_ = interpolateXY(T[celli], thermo.x(), thermo.y()); + alpha1[celli] = min(max(alpha1_, 0.0), 1.0); } alpha1.correctBoundaryConditions(); @@ -86,7 +87,7 @@ volScalarField Cp IOobject ( "Cp", - runTime.timeName(), + runTime.name(), mesh, IOobject::NO_READ, IOobject::NO_WRITE @@ -101,7 +102,7 @@ volScalarField kappa IOobject ( "kappa", - runTime.timeName(), + runTime.name(), mesh ), mesh, @@ -117,7 +118,7 @@ volScalarField p_rgh IOobject ( "p_rgh", - runTime.timeName(), + runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE @@ -132,7 +133,7 @@ volVectorField U IOobject ( "U", - runTime.timeName(), + runTime.name(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE @@ -148,22 +149,37 @@ volScalarField rhok IOobject ( "rhok", - runTime.timeName(), + runTime.name(), mesh ), 1.0 - beta*(T - Tliq) ); #include "readGravitationalAcceleration.H" -#include "readhRef.H" -#include "gh.H" + +uniformDimensionedScalarField hRef +( + IOobject + ( + "hRef", + runTime.constant(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + dimensionedScalar(dimLength, 0) +); + +dimensionedScalar ghRef(-mag(g)*hRef); +volScalarField gh("gh", (g & mesh.C()) - ghRef); +surfaceScalarField ghf("ghf", (g & mesh.Cf()) - ghRef); volScalarField p ( IOobject ( "p", - runTime.timeName(), + runTime.name(), mesh, IOobject::NO_READ, IOobject::NO_WRITE @@ -171,15 +187,12 @@ volScalarField p p_rgh + rhok*gh ); -label pRefCell = 0; -scalar pRefValue = 0.0; -setRefCell +// Create pressure reference +Foam::pressureReference pressureReference ( p, p_rgh, - pimple.dict(), - pRefCell, - pRefValue + pimple.dict() ); if (p_rgh.needReference()) @@ -188,10 +201,32 @@ if (p_rgh.needReference()) ( "p", p.dimensions(), - pRefValue - getRefCellValue(p, pRefCell) + pressureReference.refValue() + - getRefCellValue(p, pressureReference.refCell()) ); + + p_rgh = p - rhok*gh; } mesh.schemes().setFluxRequired(p_rgh.name()); #include "createMRF.H" + +autoPtr Uf; +if (mesh.dynamic() || MRF.size()) +{ + U.correctBoundaryConditions(); + + Uf = new surfaceVectorField + ( + IOobject + ( + "Uf", + runTime.name(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + fvc::interpolate(U) + ); +} diff --git a/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.C b/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.C index 8788bb7e..a0160071 100644 --- a/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.C +++ b/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.C @@ -5,7 +5,7 @@ \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2023 Oak Ridge National Laboratory + Copyright (C) 2023 Oak Ridge National Laboratory ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -91,22 +91,13 @@ Foam::marangoniFvPatchVectorField::marangoniFvPatchVectorField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::marangoniFvPatchVectorField::autoMap +void Foam::marangoniFvPatchVectorField::map ( - const fvPatchFieldMapper& m -) -{ - transformFvPatchVectorField::autoMap(m); -} - - -void Foam::marangoniFvPatchVectorField::rmap -( - const fvPatchVectorField& ptf, - const labelList& addr + const fvPatchVectorField& pvf, + const fvPatchFieldMapper& mapper ) { - transformFvPatchVectorField::rmap(ptf, addr); + transformFvPatchVectorField::map(pvf, mapper); } @@ -126,7 +117,7 @@ Foam::marangoniFvPatchVectorField::snGrad() const vectorField pif(this->patchInternalField()); // calculate the temperature gradient on the patch - const volScalarField& T = + const volScalarField& T = this->internalField().mesh().lookupObject("T"); const dimensionedScalar Tmax("Tmax", dimTemperature, Tmax_); diff --git a/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.H b/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.H index f856c091..cf9157bd 100644 --- a/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.H +++ b/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.H @@ -5,7 +5,7 @@ \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2023 Oak Ridge National Laboratory + Copyright (C) 2023 Oak Ridge National Laboratory ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -45,7 +45,7 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class marangoniFvPatchVectorField Declaration + Class marangoniFvPatchVectorField Declaration \*---------------------------------------------------------------------------*/ class marangoniFvPatchVectorField @@ -123,13 +123,12 @@ public: // Mapping functions - //- Map (and resize as needed) from self given a mapping object - // Used to update fields following mesh topology change - virtual void autoMap(const fvPatchFieldMapper&); - - //- Reverse map the given fvPatchField onto this fvPatchField - // Used to reconstruct fields - virtual void rmap(const fvPatchVectorField&, const labelList&); + //- Map the given fvPatchField onto this fvPatchField + virtual void map + ( + const fvPatchVectorField&, + const fvPatchFieldMapper& + ); //- Reset the fvPatchField to the given fvPatchField // Used for mesh to mesh mapping diff --git a/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.C b/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.C index 5f218fba..175028c6 100644 --- a/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.C +++ b/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.C @@ -2,10 +2,8 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2023 OpenFOAM Foundation \\/ M anipulation | -------------------------------------------------------------------------------- - Copyright (C) 2023 Oak Ridge National Laboratory ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -33,24 +31,6 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::mixedTemperatureFvPatchScalarField:: -mixedTemperatureFvPatchScalarField -( - const fvPatch& p, - const DimensionedField& iF -) -: - mixedFvPatchScalarField(p, iF), - h_(0.0), - emissivity_(0.0), - Tinf_(p.size(), Zero) -{ - refValue() = Zero; - refGrad() = Zero; - valueFraction() = 1; -} - - Foam::mixedTemperatureFvPatchScalarField:: mixedTemperatureFvPatchScalarField ( @@ -61,25 +41,14 @@ mixedTemperatureFvPatchScalarField : mixedFvPatchScalarField(p, iF), h_(dict.lookup("h")), - emissivity_(dict.lookup("emissivity")), + emissivity_(dict.lookupOrDefault("emissivity", 0.0)), Tinf_("Tinf", dict, p.size()) { - fvPatchScalarField::operator=(scalarField("value", dict, p.size())); + refGrad() = Zero; + valueFraction() = 0.0; - if (dict.found("refValue")) - { - // Full restart - refValue() = scalarField("refValue", dict, p.size()); - refGrad() = scalarField("refGradient", dict, p.size()); - valueFraction() = scalarField("valueFraction", dict, p.size()); - } - else - { - // Start from user entered data. Assume fixedValue. - refValue() = *this; - refGrad() = 0; - valueFraction() = 1; - } + refValue() = scalarField("Tinf", dict, p.size()); + fvPatchScalarField::operator=(refValue()); } @@ -115,28 +84,18 @@ mixedTemperatureFvPatchScalarField // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::mixedTemperatureFvPatchScalarField::autoMap -( - const fvPatchFieldMapper& m -) -{ - mixedFvPatchScalarField::autoMap(m); - m(Tinf_, Tinf_); -} - - -void Foam::mixedTemperatureFvPatchScalarField::rmap +void Foam::mixedTemperatureFvPatchScalarField::map ( const fvPatchScalarField& ptf, - const labelList& addr + const fvPatchFieldMapper& mapper ) { - mixedFvPatchScalarField::rmap(ptf, addr); + mixedFvPatchScalarField::map(ptf, mapper); const mixedTemperatureFvPatchScalarField& tiptf = refCast(ptf); - Tinf_.rmap(tiptf.Tinf_, addr); + mapper(Tinf_, tiptf.Tinf_); } @@ -161,7 +120,7 @@ void Foam::mixedTemperatureFvPatchScalarField::updateCoeffs() return; } - refValue() = Tinf_; + mixedFvPatchScalarField::refValue() = (Tinf_); const scalarField& Tp(*this); @@ -169,18 +128,14 @@ void Foam::mixedTemperatureFvPatchScalarField::updateCoeffs() patch().lookupPatchField("kappa"); const scalar sigma_(5.67e-8); - + scalarField hEff_ ( - h_ + sigma_*emissivity_*(sqr(Tp) + sqr(Tinf_))*(Tp + Tinf_) + h_ + sigma_ * emissivity_ * (sqr(Tp) + sqr(Tinf_)) * (Tp + Tinf_) ); - valueFraction() = - 1.0/ - ( - 1.0 - + kappa_*patch().deltaCoeffs()/max(hEff_, 1e-15) - ); + valueFraction() = + 1.0 / (1.0 + kappa_ * patch().deltaCoeffs() / max(hEff_, 1e-15)); refGrad() = Zero; diff --git a/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.H b/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.H index 3d813a34..68ed7a61 100644 --- a/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.H +++ b/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.H @@ -2,10 +2,8 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2023 OpenFOAM Foundation \\/ M anipulation | -------------------------------------------------------------------------------- - Copyright (C) 2023 Oak Ridge National Laboratory ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -30,15 +28,40 @@ Description This boundary condition provides a mixedTemperature condition, calculated as: -SourceFiles - mixedTemperatureFvPatchScalarField.C + \f[ + Q = h * (T - T_inf) + e * sigma * (T^4 - T_inf^$) + \f] + + where + \vartable + h | Convection coefficient [W/m^2/K] + e | Surface emissivity [-] + sigma | Stefan-Boltzmann constant [5.67e-8 W/m^2/K^4] + \endtable + +Usage + \table + Property | Description | Req'd? | Default + h | Uniform convection coeff | yes | + emissivity | Uniform surface emissivity | no | 0.0 + Tinf | Ambient temperature | yes | + \endtable + + Example of the boundary condition specification: + \verbatim + + { + type mixedTemperature; + h 10.0; + emissivity 0.4; + } + \endverbatim SourceFiles mixedTemperatureFvPatchScalarField.C \*---------------------------------------------------------------------------*/ - #ifndef BC_H #define BC_H @@ -60,13 +83,13 @@ class mixedTemperatureFvPatchScalarField { // Private Data - //- convective heat transfer coefficient [W/m^2/k] + //- Convective heat transfer coefficient [W/m^2/K] scalar h_; - //- effective emissivity of the boundary + //- Effective emissivity of the boundary scalar emissivity_; - //- ambient temperature + //- Ambient temperature scalarField Tinf_; @@ -78,13 +101,6 @@ public: // Constructors - //- Construct from patch and internal field - mixedTemperatureFvPatchScalarField - ( - const fvPatch&, - const DimensionedField& - ); - //- Construct from patch, internal field and dictionary mixedTemperatureFvPatchScalarField ( @@ -137,13 +153,12 @@ public: // Mapping functions - //- Map (and resize as needed) from self given a mapping object - // Used to update fields following mesh topology change - virtual void autoMap(const fvPatchFieldMapper&); - - //- Reverse map the given fvPatchField onto this fvPatchField - // Used to reconstruct fields - virtual void rmap(const fvPatchScalarField&, const labelList&); + //- Map the given fvPatchField onto this fvPatchField + virtual void map + ( + const fvPatchScalarField&, + const fvPatchFieldMapper& + ); //- Reset the fvPatchField to the given fvPatchField // Used for mesh to mesh mapping diff --git a/applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.C b/applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.C index f8878348..556a7c4c 100644 --- a/applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.C +++ b/applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.C @@ -5,7 +5,7 @@ \\ / A nd | Copyright (C) 2024 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2023 Oak Ridge National Laboratory + Copyright (C) 2023 Oak Ridge National Laboratory ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -29,20 +29,26 @@ License #include "Time.H" #include "fvMesh.H" #include "addToRunTimeSelectionTable.H" -#include "volFields.H" +#include "volPointInterpolation.H" #include "OFstream.H" #include "OSspecific.H" -#include "labelVector.H" #include "pointMVCWeight.H" +#include "Pstream.H" -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // +#include "polyTopoChangeMap.H" +#include "polyMeshMap.H" +#include "polyDistributionMap.H" + +#include + +// * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * * * // namespace Foam { namespace functionObjects { defineTypeNameAndDebug(ExaCA, 0); - + addToRunTimeSelectionTable ( functionObject, @@ -52,6 +58,7 @@ namespace functionObjects } } + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::functionObjects::ExaCA::ExaCA @@ -62,25 +69,16 @@ Foam::functionObjects::ExaCA::ExaCA ) : fvMeshFunctionObject(name, runTime, dict), - T_(mesh_.lookupObject>("T")), - vpi_(mesh_), - Tp_ - ( - IOobject - ( - "Tp_", - runTime.timeName(), - mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE - ), - vpi_.interpolate(T_) - ), - searchEngine_(mesh_, polyMesh::CELL_TETS) + T_(mesh_.lookupObject("T")), + meshChanged_(true), + nPoints_(labelVector(vector::one)), + box_(point::max, point::min), + isoValue_(Zero), + dx_(Zero) { read(dict); - - setOverlapCells(); + + updateMeshState(); } @@ -90,22 +88,332 @@ Foam::functionObjects::ExaCA::~ExaCA() {} +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + +void Foam::functionObjects::ExaCA::resetMeshState() +{ + cellMapIds.clear(); + + overlapCells.clear(); + + meshChanged_ = true; +} + + +void Foam::functionObjects::ExaCA::updateMeshState() +{ + cellMapIds.clear(); + + overlapCells.clear(); + + setOverlapCells(); + + meshChanged_ = false; +} + + +Foam::label Foam::functionObjects::ExaCA::pointId +( + const label i, + const label j, + const label k +) const +{ + return i + nPoints_.x()*(j + nPoints_.y()*k); +} + + +Foam::label Foam::functionObjects::ExaCA::getCellMap(const label celli) +{ + if (cellMapIds.found(celli)) + { + return cellMapIds[celli]; + } + + const label mapId = cellMaps.size(); + + cellMapIds.set(celli, mapId); + + cellMaps.append(makeCellMap(celli)); + + return mapId; +} + + +Foam::functionObjects::ExaCA::cellMap +Foam::functionObjects::ExaCA::makeCellMap(const label celli) const +{ + cellMap cm; + + const pointField& points = mesh_.points(); + + const labelList& vertices = mesh_.cellPoints()[celli]; + + boundBox cellBb(point::max, point::min); + + const vector extend = 1e-10*vector::one; + + forAll(vertices, vertexi) + { + const point& p = points[vertices[vertexi]]; + + cellBb.min() = min(cellBb.min(), p - extend); + cellBb.max() = max(cellBb.max(), p + extend); + } + + if (!cellBb.overlaps(box_)) + { + return cm; + } + + const point& boxMax = box_.max(); + + const label i0 = + max + ( + label(0), + label(std::ceil((boxMax.x() - cellBb.max().x())/dx_)) + ); + + const label i1 = + min + ( + nPoints_.x() - 1, + label(std::floor((boxMax.x() - cellBb.min().x())/dx_)) + ); + + const label j0 = + max + ( + label(0), + label(std::ceil((boxMax.y() - cellBb.max().y())/dx_)) + ); + + const label j1 = + min + ( + nPoints_.y() - 1, + label(std::floor((boxMax.y() - cellBb.min().y())/dx_)) + ); + + const label k0 = + max + ( + label(0), + label(std::ceil((boxMax.z() - cellBb.max().z())/dx_)) + ); + + const label k1 = + min + ( + nPoints_.z() - 1, + label(std::floor((boxMax.z() - cellBb.min().z())/dx_)) + ); + + if (i1 < i0 || j1 < j0 || k1 < k0) + { + return cm; + } + + const cell& cFaces = mesh_.cells()[celli]; + + const label nPlanes = cFaces.size(); + + scalarField nx(nPlanes); + scalarField ny(nPlanes); + scalarField nz(nPlanes); + scalarField d(nPlanes); + + const vectorField& faceAreas = mesh_.faceAreas(); + const pointField& faceCentres = mesh_.faceCentres(); + const labelUList& faceOwner = mesh_.faceOwner(); + + forAll(cFaces, cFacei) + { + const label facei = cFaces[cFacei]; + + const vector& Sf = faceAreas[facei]; + const point& Cf = faceCentres[facei]; + + if (faceOwner[facei] == celli) + { + nx[cFacei] = Sf.x(); + ny[cFacei] = Sf.y(); + nz[cFacei] = Sf.z(); + } + else + { + nx[cFacei] = -Sf.x(); + ny[cFacei] = -Sf.y(); + nz[cFacei] = -Sf.z(); + } + + d[cFacei] = + nx[cFacei]*Cf.x() + + ny[cFacei]*Cf.y() + + nz[cFacei]*Cf.z(); + } + + const scalar xMax = boxMax.x(); + const scalar yMax = boxMax.y(); + const scalar zMax = boxMax.z(); + + const scalar eps = 1e-10; + + for (label k = k0; k <= k1; ++k) + { + const scalar z = zMax - scalar(k)*dx_; + const scalar sz = z - eps; + + for (label j = j0; j <= j1; ++j) + { + const scalar y = yMax - scalar(j)*dx_; + const scalar sy = y - eps; + + for (label i = i0; i <= i1; ++i) + { + const scalar x = xMax - scalar(i)*dx_; + const scalar sx = x - eps; + + bool pointInCell = true; + + for (label planei = 0; planei < nPlanes; ++planei) + { + if + ( + nx[planei]*sx + + ny[planei]*sy + + nz[planei]*sz + - d[planei] + > 0 + ) + { + pointInCell = false; + break; + } + } + + if (!pointInCell) + { + continue; + } + + const point pt(x, y, z); + + const point spt(sx, sy, sz); + + pointMVCWeight cpw(mesh_, spt, celli); + + cm.pointIds.append(pointId(i, j, k)); + + cm.positions.append(pt); + + cm.weights.append(cpw.weights()); + } + } + } + + cm.pointIds.shrink(); + + cm.positions.shrink(); + + cm.weights.shrink(); + + return cm; +} + + +void Foam::functionObjects::ExaCA::appendInterval +( + const label celli, + const pointScalarField& Tp0, + const pointScalarField& Tp, + const scalar t0, + const scalar t +) +{ + const label mapId = getCellMap(celli); + + if (!cellMaps[mapId].pointIds.size()) + { + return; + } + + const labelList& vertices = mesh_.cellPoints()[celli]; + + eventInterval event; + + event.mapId = mapId; + + event.t0 = t0; + + event.t1 = t; + + event.psi0.setSize(vertices.size()); + + event.psi.setSize(vertices.size()); + + forAll(vertices, pI) + { + event.psi0[pI] = Tp0[vertices[pI]]; + + event.psi[pI] = Tp[vertices[pI]]; + } + + intervals.append(event); +} + + +void Foam::functionObjects::ExaCA::writeData() const +{ + const fileName exacaPath + ( + mesh_.time().rootPath()/mesh_.time().globalCaseName()/"ExaCA" + ); + + mkDir(exacaPath); + + OFstream os + ( + exacaPath + "/" + "data_" + Foam::name(Pstream::myProcNo()) + ".csv" + ); + + os << "x,y,z,tm,ts,cr" << endl; + + forAll(data, i) + { + const List& row = data[i]; + + os << row[0] << "," + << row[1] << "," + << row[2] << "," + << row[3] << "," + << row[4] << "," + << row[5] << "\n"; + } +} + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // bool Foam::functionObjects::ExaCA::read(const dictionary& dict) { box_ = dict.lookup("box"); - + isoValue_ = dict.lookup("isoValue"); - - dx_ = dict.lookup("dx"); + + dx_ = dict.lookup("dx"); + + nPoints_ = labelVector(vector::one + box_.span()/dx_); return true; } + void Foam::functionObjects::ExaCA::setOverlapCells() { - // create a compact cell-stencil using the overlap sub-space + overlapCells.clear(); + const pointField& points = mesh_.points(); boundBox procBb(points); @@ -139,237 +447,160 @@ void Foam::functionObjects::ExaCA::setOverlapCells() Foam::wordList Foam::functionObjects::ExaCA::fields() const { - return wordList::null(); + return wordList(1, T_.name()); } bool Foam::functionObjects::ExaCA::execute() -{ - const pointScalarField Tp0_("Tp0_", Tp_); +{ + if (meshChanged_) + { + updateMeshState(); + } + + const volPointInterpolation& vpi = volPointInterpolation::New(mesh_); + + const pointScalarField Tp0 + ( + IOobject + ( + "Tp0_", + mesh_.time().name(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + vpi.interpolate(T_.oldTime()) + ); + + const pointScalarField Tp + ( + IOobject + ( + "Tp_", + mesh_.time().name(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + vpi.interpolate(T_) + ); + + const scalar t = mesh_.time().value(); - Tp_ = vpi_.interpolate(T_); + const scalar t0 = t - mesh_.time().deltaTValue(); - // capture events at interface cells: {cell id, time, vertex temperatures} - const scalar t_ = mesh_.time().value(); - const scalar t0_ = t_ - mesh_.time().deltaTValue(); - forAll(overlapCells, i) { - label celli = overlapCells[i]; + const label celli = overlapCells[i]; const labelList& vertices = mesh_.cellPoints()[celli]; label c0 = 0; + label c1 = 0; forAll(vertices, pI) { - if (Tp0_[vertices[pI]] >= isoValue_) + if (Tp0[vertices[pI]] >= isoValue_) { c0++; } - if (Tp_[vertices[pI]] >= isoValue_) + if (Tp[vertices[pI]] >= isoValue_) { c1++; } } - + const label n = vertices.size(); - // overshoot correction: interface jumped this cell during time step - if ( !(c0 % n) && !(c1 % n) ) + if (!(c0 % n) && !(c1 % n)) { - if (c0 != c1) + if (c0 != c1) { c0 = 0; + c1 = 1; } } - - // capture solidification events + c0 %= n; + c1 %= n; - + if (c0 || c1) { - List event(n + 2); - - event[0] = celli; - - // add previous event - if (c1 && !c0) - { - event[1] = t0_; - - forAll(vertices, pI) - { - event[pI + 2] = Tp0_[vertices[pI]]; - } - - events.append(event); - } - - // add current event - event[1] = t_; - - forAll(vertices, pI) - { - event[pI + 2] = Tp_[vertices[pI]]; - } - - events.append(event); + appendInterval(celli, Tp0, Tp, t0, t); } } - + return true; } -void Foam::functionObjects::ExaCA::mapPoints() + +void Foam::functionObjects::ExaCA::interpolate() { - if (events.size() == 0) + if (intervals.size() == 0) { + writeData(); + return; } - // find event sub-space before constructing interpolant weights - const pointField& points = mesh_.points(); - - const vector extend = 1e-10*vector::one; + meltTimes.clear(); - boundBox eventBb(point::max, point::min); + data.clear(); - for (label i = 1; i < events.size(); i++) + forAll(intervals, eventi) { - const label celli = events[i][0]; + const eventInterval& event = intervals[eventi]; - if (celli == events[i - 1][0]) + if (event.t1 <= event.t0 + small) { continue; } - const labelList& vertices = mesh_.cellPoints()[celli]; + const cellMap& cm = cellMaps[event.mapId]; - forAll(vertices, i) + forAll(cm.pointIds, pointi) { - eventBb.min() = min(eventBb.min(), points[vertices[i]] - extend); - eventBb.max() = max(eventBb.max(), points[vertices[i]] + extend); - } - } - - // find interpolant weigthts for each point - label pI = 0; - label seedi = events[0][0]; - - pointsInCell.setSize(mesh_.nCells()); + scalar tp0 = Zero; - const labelVector nPoints(vector::one + box_.span() / dx_); + scalar tp = Zero; - Info << "starting point loop" << endl; + const scalarField& w = cm.weights[pointi]; - for (label k=0; k < nPoints.z(); ++k) - { - for (label j=0; j < nPoints.y(); ++j) - { - for (label i=0; i < nPoints.x(); ++i) + forAll(w, j) { - const point pt = box_.max() - vector(i, j, k)*dx_; - - if (eventBb.contains(pt)) - { - // shift point during search to avoid edges in pointMVC - const point spt = pt - vector::one*1e-10; + tp0 += w[j]*event.psi0[j]; - label celli = searchEngine_.findCell(spt, seedi, true); - - if (celli != -1) - { - positions.append(pt); - - pointMVCWeight cpw(mesh_, spt, celli); - - weights.append(cpw.weights()); - - pointsInCell[celli].append(pI); - - pI++; - } - - seedi = celli; - } + tp += w[j]*event.psi[j]; } - } - } - - positions.shrink(); - - weights.shrink(); - - for (auto& pic : pointsInCell) - { - pic.shrink(); - } -} - -void Foam::functionObjects::ExaCA::interpolate() -{ - if (events.size() == 0) - { - return; - } - - // format events in ExaCA reduced data format - DynamicList> data; - - List tm; - tm.setSize(pointsInCell[events[0][0]].size(), events[0][1]); - for (label i = 1; i < events.size(); i++) - { - const List prevEvent = events[i - 1]; - - const List currEvent = events[i]; - - label celli = currEvent[0]; - - if (currEvent[0] != prevEvent[0]) - { - tm.setSize(pointsInCell[celli].size(), currEvent[1]); - continue; - } - - scalar prevTime = prevEvent[1]; - - scalar currTime = currEvent[1]; - - List psi0(prevEvent.size() - 2); - List psi(currEvent.size() - 2); - - for (int j = 0; j < psi.size(); j++) - { - psi0[j] = prevEvent[j + 2]; - psi[j] = currEvent[j + 2]; - } - - int p = 0; - for (const label& pointi : pointsInCell[celli]) - { - scalar tp0 = Zero; - scalar tp = Zero; - - List w = weights[pointi]; - - for (int j = 0; j < psi.size(); j++) - { - tp0 += w[j]*psi0[j]; - tp += w[j]*psi[j]; - } + const label caPointId = cm.pointIds[pointi]; if ((tp <= isoValue_) && (tp0 > isoValue_)) { - const point& pt = positions[pointi]; - - scalar m_ = min(max((isoValue_ - tp0)/(tp - tp0), 0), 1); + const point& pt = cm.positions[pointi]; + + const scalar m = + min + ( + max + ( + (isoValue_ - tp0)/(tp - tp0), + scalar(0) + ), + scalar(1) + ); + + scalar tm = event.t0; + + if (meltTimes.found(caPointId)) + { + tm = meltTimes[caPointId]; + } data.append ( @@ -377,80 +608,58 @@ void Foam::functionObjects::ExaCA::interpolate() pt[0], pt[1], pt[2], - tm[p], - prevTime + m_*(currTime - prevTime), - (tp0 - tp) / (currTime - prevTime) + tm, + event.t0 + m*(event.t1 - event.t0), + (tp0 - tp)/(event.t1 - event.t0) } ); } else if ((tp > isoValue_) && (tp0 <= isoValue_)) { - scalar m_ = min(max((isoValue_ - tp0)/(tp - tp0), 0), 1); - - tm[p] = prevTime + m_*(currTime - prevTime); + const scalar m = + min + ( + max + ( + (isoValue_ - tp0)/(tp - tp0), + scalar(0) + ), + scalar(1) + ); + + meltTimes.set(caPointId, event.t0 + m*(event.t1 - event.t0)); } - - p++; } } - // write the events for each processor to their own file - const fileName exacaPath - ( - mesh_.time().rootPath()/mesh_.time().globalCaseName()/"ExaCA" - ); + data.shrink(); - OFstream os - ( - exacaPath + "/" + "data_" + Foam::name(Pstream::myProcNo()) + ".csv" - ); - - os << "x,y,z,tm,ts,cr" << endl; - - for(int i=0; i < data.size(); i++) - { - int n = data[i].size()-1; - - for(int j=0; j < n; j++) - { - os << data[i][j] << ","; - } - os << data[i][n] << "\n"; - } + writeData(); } + bool Foam::functionObjects::ExaCA::end() { - //- sort events by cell and in time - events.shrink(); + intervals.shrink(); - sort(events); + cellMaps.shrink(); - Info<< "Number of solidification events: " - << returnReduce(events.size(), sumOp()) << endl; + Info<< "Number of ExaCA event intervals: " + << returnReduce(intervals.size(), sumOp