diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 626c2dd..bd25471 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -19,7 +19,7 @@ jobs:
         cxx: ['g++']
         cmake_build_type: ['Release']
        
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-latest
     container:
       image: docker.io/openfoam/openfoam10-paraview510
       options: --user root
@@ -29,11 +29,8 @@ jobs:
       - name: Build AdditiveFOAM
         run: |
           . /opt/openfoam10/etc/bashrc
-          cd applications/solvers/additiveFoam/movingHeatSource
-          wmake libso
-          cd -
           cd applications/solvers/additiveFoam
-          wmake
+          ./Allwmake
       - name: Test AdditiveFOAM
         run: |
           . /opt/openfoam10/etc/bashrc
diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml
new file mode 100644
index 0000000..32afd36
--- /dev/null
+++ b/.github/workflows/draft-pdf.yml
@@ -0,0 +1,24 @@
+name: Draft PDF
+on: [push]
+
+jobs:
+  paper:
+    runs-on: ubuntu-latest
+    name: Paper Draft
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+      - name: Build draft PDF
+        uses: openjournals/openjournals-draft-action@master
+        with:
+          journal: joss
+          # This should be the path to the paper within your repo.
+          paper-path: doc/paper.md
+      - name: Upload
+        uses: actions/upload-artifact@v4
+        with:
+          name: paper
+          # This is the output path where Pandoc will write the compiled
+          # PDF. Note, this should be the same directory as the input
+          # paper.md
+          path: doc/paper.pdf
diff --git a/README.md b/README.md
index 7cbc323..4b911d3 100644
--- a/README.md
+++ b/README.md
@@ -1,5 +1,12 @@
-# AdditiveFOAM
-AdditiveFOAM is a free, open source heat and mass transfer software for simulations of Additive Manufacturing (AM) released by Oak Ridge National Laboratory. It is built upon OpenFOAM, a free, open source computational fluid dynamics (CFD) software package released by the OpenFOAM Foundation.
+<img src="https://raw.githubusercontent.com/ORNL-MDF/additivefoam-website-assets/main/images/logo.png" alt="image" width="800">
+
+---
+`AdditiveFOAM` is a free, open source heat and mass transfer software for simulations of Additive Manufacturing (AM) released by Oak Ridge National Laboratory. It is built upon OpenFOAM, a free, open source computational fluid dynamics (CFD) software package released by the OpenFOAM Foundation.
+
+## Documentation
+[![Documentation Status][docs-badge]][docs-url]
+
+The documentation for `AdditiveFOAM` is hosted on [GitHub Pages](https://ornl.github.io/AdditiveFOAM/).
 
 ### Repository Features
 | Link                                                | Description                              |
@@ -8,11 +15,6 @@ AdditiveFOAM is a free, open source heat and mass transfer software for simulati
 | [utilities](applications/utilities)                       | Utilities for post-processing and code wrappers |
 | [tutorials](tutorials)                                     | Tutorial cases |
 
-## Documentation
-[![Documentation Status][docs-badge]][docs-url]
-
-The documentation for `AdditiveFOAM` is hosted on [GitHub Pages](https://ornl.github.io/AdditiveFOAM/).
-
 ## Installation and Dependencies
 [![OpenFOAM-10](https://img.shields.io/badge/OpenFOAM-10-blue.svg)](https://github.com/OpenFOAM/OpenFOAM-10)
 
diff --git a/applications/solvers/additiveFoam/Allwclean b/applications/solvers/additiveFoam/Allwclean
new file mode 100755
index 0000000..20a96ae
--- /dev/null
+++ b/applications/solvers/additiveFoam/Allwclean
@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+wclean libso functionObjects/ExaCA
+wclean libso movingHeatSource
+wclean
+
+#------------------------------------------------------------------------------
diff --git a/applications/solvers/additiveFoam/Allwmake b/applications/solvers/additiveFoam/Allwmake
new file mode 100755
index 0000000..2e2f39d
--- /dev/null
+++ b/applications/solvers/additiveFoam/Allwmake
@@ -0,0 +1,11 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Parse arguments for library compilation
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+
+wmake $targetType functionObjects/ExaCA
+wmake $targetType movingHeatSource
+wmake $targetType
+
+#------------------------------------------------------------------------------
diff --git a/applications/solvers/additiveFoam/Make/files b/applications/solvers/additiveFoam/Make/files
index 5611e03..10e328d 100644
--- a/applications/solvers/additiveFoam/Make/files
+++ b/applications/solvers/additiveFoam/Make/files
@@ -1,5 +1,5 @@
-BCs/mixedTemperature/mixedTemperatureFvPatchScalarField.C
-BCs/marangoni/marangoniFvPatchVectorField.C
+derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.C
+derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.C
 
 additiveFoam.C
 
diff --git a/applications/solvers/additiveFoam/additiveFoam.C b/applications/solvers/additiveFoam/additiveFoam.C
index 98d64ca..b5769e7 100644
--- a/applications/solvers/additiveFoam/additiveFoam.C
+++ b/applications/solvers/additiveFoam/additiveFoam.C
@@ -38,10 +38,9 @@ Description
 #include "Polynomial.H"
 
 #include "interpolateXY/interpolateXY.H"
-
 #include "movingHeatSourceModel.H"
-#include "foamToExaCA/foamToExaCA.H"
 
+#include "EulerDdtScheme.H"
 #include "CrankNicolsonDdtScheme.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -65,8 +64,6 @@ int main(int argc, char *argv[])
 
     scalar alphaCoNum = 0.0;
 
-    foamToExaCA ExaCA(T);
-
     movingHeatSourceModel sources(mesh);
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -97,8 +94,6 @@ int main(int argc, char *argv[])
 
         #include "thermo/TEqn.H"
         
-        ExaCA.update();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
@@ -106,8 +101,6 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    ExaCA.write();
-
     return 0;
 }
 
diff --git a/applications/solvers/additiveFoam/createFields.H b/applications/solvers/additiveFoam/createFields.H
index 2b14874..7e159dc 100644
--- a/applications/solvers/additiveFoam/createFields.H
+++ b/applications/solvers/additiveFoam/createFields.H
@@ -21,7 +21,7 @@ volScalarField alpha3
         runTime.timeName(),
         mesh,
         IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
+        IOobject::AUTO_WRITE
     ),
     mesh,
     dimensionedScalar("alpha.powder", dimless, 0.0),
diff --git a/applications/solvers/additiveFoam/BCs/marangoni/marangoniFvPatchVectorField.C b/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.C
similarity index 100%
rename from applications/solvers/additiveFoam/BCs/marangoni/marangoniFvPatchVectorField.C
rename to applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.C
diff --git a/applications/solvers/additiveFoam/BCs/marangoni/marangoniFvPatchVectorField.H b/applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.H
similarity index 100%
rename from applications/solvers/additiveFoam/BCs/marangoni/marangoniFvPatchVectorField.H
rename to applications/solvers/additiveFoam/derivedFvPatchFields/marangoni/marangoniFvPatchVectorField.H
diff --git a/applications/solvers/additiveFoam/BCs/mixedTemperature/mixedTemperatureFvPatchScalarField.C b/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.C
similarity index 100%
rename from applications/solvers/additiveFoam/BCs/mixedTemperature/mixedTemperatureFvPatchScalarField.C
rename to applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.C
diff --git a/applications/solvers/additiveFoam/BCs/mixedTemperature/mixedTemperatureFvPatchScalarField.H b/applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.H
similarity index 100%
rename from applications/solvers/additiveFoam/BCs/mixedTemperature/mixedTemperatureFvPatchScalarField.H
rename to applications/solvers/additiveFoam/derivedFvPatchFields/mixedTemperature/mixedTemperatureFvPatchScalarField.H
diff --git a/applications/solvers/additiveFoam/foamToExaCA/foamToExaCA.C b/applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.C
similarity index 65%
rename from applications/solvers/additiveFoam/foamToExaCA/foamToExaCA.C
rename to applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.C
index d7ae22e..d1cb1ed 100644
--- a/applications/solvers/additiveFoam/foamToExaCA/foamToExaCA.C
+++ b/applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  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) 2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                 Copyright (C) 2023 Oak Ridge National Laboratory                
@@ -25,10 +25,13 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "foamToExaCA/foamToExaCA.H"
+#include "ExaCA.H"
+#include "Time.H"
+#include "fvMesh.H"
 #include "addToRunTimeSelectionTable.H"
-#include "DynamicList.H"
-#include "interpolation.H"
+#include "volFields.H"
+#include "OFstream.H"
+#include "OSspecific.H"
 #include "labelVector.H"
 #include "pointMVCWeight.H"
 
@@ -36,86 +39,84 @@ License
 
 namespace Foam
 {
-    defineTypeNameAndDebug(foamToExaCA, 0);
+namespace functionObjects
+{
+    defineTypeNameAndDebug(ExaCA, 0);
+    
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        ExaCA,
+        dictionary
+    );
+}
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::foamToExaCA::foamToExaCA
+Foam::functionObjects::ExaCA::ExaCA
 (
-    const volScalarField& T
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
 )
 :
-    IOdictionary
-    (
-        IOobject
-        (
-            "foamToExaCADict",
-            T.time().constant(),
-            T.mesh(),
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        )
-    ),
-
-    mesh_(T.mesh()),
-
-    runTime_(T.time()),
-
-    T_(T),
-
+    fvMeshFunctionObject(name, runTime, dict),
+    T_(mesh_.lookupObject<VolField<scalar>>("T")),
     vpi_(mesh_),
-
     Tp_
     (
         IOobject
         (
             "Tp_",
-            runTime_.timeName(),
+            runTime.timeName(),
             mesh_,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         vpi_.interpolate(T_)
     ),
-
-    execute_(false)
+    searchEngine_(mesh_, polyMesh::CELL_TETS)
 {
-    if (this->headerOk())
-    {
-        execute_ = this->lookup<Switch>("execute");
-    }
-
-    initialize();
+    read(dict);
+    
+    setOverlapCells();
 }
 
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::foamToExaCA::initialize()
-{
-    if (!execute_)
-    {
-        return;
-    }
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::ExaCA::~ExaCA()
+{}
+
 
-    box_ = this->lookup("box");
-    iso_ = this->lookup<scalar>("isotherm");
-    dx_  = this->lookup<scalar>("dx");
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-    treeBoundBox bb(box_);
+bool Foam::functionObjects::ExaCA::read(const dictionary& dict)
+{
+    box_ = dict.lookup("box");
+    
+    isoValue_ = dict.lookup<scalar>("isoValue");
+    
+    dx_  = dict.lookup<scalar>("dx");
+
+    return true;
+}
 
+void Foam::functionObjects::ExaCA::setOverlapCells()
+{
     // create a compact cell-stencil using the overlap sub-space
     const pointField& points = mesh_.points();
 
-    treeBoundBox procBb(points);
+    boundBox procBb(points);
 
     const vector extend = 1e-10*vector::one;
 
-    if (procBb.overlaps(bb))
+    if (procBb.overlaps(box_))
     {
         forAll(mesh_.cells(), celli)
         {
-            treeBoundBox cellBb(point::max, point::min);
+            boundBox cellBb(point::max, point::min);
 
             const labelList& vertices = mesh_.cellPoints()[celli];
 
@@ -125,78 +126,81 @@ void Foam::foamToExaCA::initialize()
                 cellBb.max() = max(cellBb.max(), points[vertices[i]] + extend);
             }
 
-            if (cellBb.overlaps(bb))
+            if (cellBb.overlaps(box_))
             {
-                compactCells.append(celli);
+                overlapCells.append(celli);
             }
         }
     }
 
-    compactCells.shrink();
+    overlapCells.shrink();
 }
 
-void Foam::foamToExaCA::update()
+
+Foam::wordList Foam::functionObjects::ExaCA::fields() const
 {
-    if (!execute_)
-    {
-        return;
-    }
+    return wordList::null();
+}
+
 
+bool Foam::functionObjects::ExaCA::execute()
+{    
     const pointScalarField Tp0_("Tp0_", Tp_);
 
     Tp_ = vpi_.interpolate(T_);
 
-    // capture events at interface cells: {cell id, time, vertice temperatures}
-    forAll(compactCells, i)
+    // 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 = compactCells[i];
+        label celli = overlapCells[i];
 
         const labelList& vertices = mesh_.cellPoints()[celli];
 
-        label prev = 0;
-        label curr = 0;
+        label c0 = 0;
+        label c1 = 0;
 
         forAll(vertices, pI)
         {
-            if (Tp0_[vertices[pI]] >= iso_)
+            if (Tp0_[vertices[pI]] >= isoValue_)
             {
-                prev++;
+                c0++;
             }
 
-            if (Tp_[vertices[pI]] >= iso_)
+            if (Tp_[vertices[pI]] >= isoValue_)
             {
-                curr++;
+                c1++;
             }
         }
-
-        prev = prev % vertices.size();
-        curr = curr % vertices.size();
+        
+        const label n = vertices.size();
 
         // overshoot correction: interface jumped this cell during time step
-        if (!prev && !curr)
+        if ( !(c0 % n) && !(c1 % n) )
         {
-            if
-            (
-                (T_[celli] > iso_ && T_.oldTime()[celli] < iso_)
-             || (T_[celli] < iso_ && T_.oldTime()[celli] > iso_)
-            )
+            if  (c0 != c1)
             {
-                prev = 0;
-                curr = 1;
+                c0 = 0;
+                c1 = 1;
             }
         }
-
-        // capture the solidification events
-        if (prev || curr)
+        
+        // capture solidification events
+        c0 %= n;
+        c1 %= n;
+        
+        if (c0 || c1)
         {
-            List<scalar> event(vertices.size() + 2);
+            List<scalar> event(n + 2);
 
             event[0] = celli;
         
             // add previous event
-            if (curr && !prev)
+            if (c1 && !c0)
             {
-                event[1] = runTime_.value() - runTime_.deltaTValue();
+                event[1] = t0_;
 
                 forAll(vertices, pI)
                 {
@@ -207,7 +211,7 @@ void Foam::foamToExaCA::update()
             }
 
             // add current event
-            event[1] = runTime_.value();
+            event[1] = t_;
 
             forAll(vertices, pI)
             {
@@ -217,25 +221,115 @@ void Foam::foamToExaCA::update()
             events.append(event);
         }
     }
+    
+    return true;
 }
 
-void Foam::foamToExaCA::interpolatePoints()
+void Foam::functionObjects::ExaCA::mapPoints()
 {
-    // dynamic list for exaca reduced data format
+    if (events.size() == 0)
+    {
+        return;
+    }
+
+    // find event sub-space before constructing interpolant weights
+    const pointField& points = mesh_.points();
+
+    const vector extend = 1e-10*vector::one;
+
+    boundBox eventBb(point::max, point::min);
+
+    for (label i = 1; i < events.size(); i++)
+    {
+        const label celli = events[i][0];
+
+        if (celli == events[i - 1][0])
+        {
+            continue;
+        }
+
+        const labelList& vertices = mesh_.cellPoints()[celli];
+
+        forAll(vertices, i)
+        {
+            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());
+
+    const labelVector nPoints(vector::one + box_.span() / dx_);
+
+    Info << "starting point loop" << endl;
+
+    for (label k=0; k < nPoints.z(); ++k)
+    {
+        for (label j=0; j < nPoints.y(); ++j)
+        {
+            for (label i=0; i < nPoints.x(); ++i)
+            {
+                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;
+
+                    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;
+                }
+            }
+        }
+    }
+
+    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<List<scalar>> data;
 
-    // initialize melting time for first event
     List<scalar> tm;
     tm.setSize(pointsInCell[events[0][0]].size(), events[0][1]);
 
-    // events are order by cell id, and each cell id is ordered in time
     for (label i = 1; i < events.size(); i++)
     {
         const List<scalar> prevEvent = events[i - 1];
 
         const List<scalar> currEvent = events[i];
 
-        // first event for the cell. set the melting time and continue
         label celli = currEvent[0];
 
         if (currEvent[0] != prevEvent[0])
@@ -244,7 +338,6 @@ void Foam::foamToExaCA::interpolatePoints()
             continue;
         }
 
-        // extract event information
         scalar prevTime = prevEvent[1];
 
         scalar currTime = currEvent[1];
@@ -272,11 +365,11 @@ void Foam::foamToExaCA::interpolatePoints()
                 tp  += w[j]*psi[j];
             }
 
-            if ((tp <= iso_) && (tp0 > iso_))
+            if ((tp <= isoValue_) && (tp0 > isoValue_))
             {
                 const point& pt = positions[pointi];
 
-                scalar m_ = min(max((iso_ - tp0)/(tp - tp0), 0), 1);
+                scalar m_ = min(max((isoValue_ - tp0)/(tp - tp0), 0), 1);
 
                 data.append
                 (
@@ -290,9 +383,9 @@ void Foam::foamToExaCA::interpolatePoints()
                     }
                 );
             }
-            else if ((tp > iso_) && (tp0 <= iso_))
+            else if ((tp > isoValue_) && (tp0 <= isoValue_))
             {
-                scalar m_ = min(max((iso_ - tp0)/(tp - tp0), 0), 1);
+                scalar m_ = min(max((isoValue_ - tp0)/(tp - tp0), 0), 1);
 
                 tm[p] = prevTime + m_*(currTime - prevTime);
             }
@@ -304,7 +397,7 @@ void Foam::foamToExaCA::interpolatePoints()
     // write the events for each processor to their own file
     const fileName exacaPath
     (
-        runTime_.rootPath()/runTime_.globalCaseName()/"ExaCA"
+        mesh_.time().rootPath()/mesh_.time().globalCaseName()/"ExaCA"
     );
 
     OFstream os
@@ -324,122 +417,47 @@ void Foam::foamToExaCA::interpolatePoints()
     }
 }
 
-
-void Foam::foamToExaCA::mapPoints(const meshSearch& searchEngine)
+bool Foam::functionObjects::ExaCA::end()
 {
-    // find event sub-space before constructing interpolants
-    const pointField& points = mesh_.points();
-
-    const vector extend = 1e-10*vector::one;
-
-    treeBoundBox eventBb(point::max, point::min);
-
-    for (label i = 1; i < events.size(); i++)
-    {
-        const label celli = events[i][0];
-
-        if (celli == events[i - 1][0])
-        {
-            continue;
-        }
-
-        const labelList& vertices = mesh_.cellPoints()[celli];
-
-        forAll(vertices, i)
-        {
-            eventBb.min() = min(eventBb.min(), points[vertices[i]] - extend);
-            eventBb.max() = max(eventBb.max(), points[vertices[i]] + extend);
-        }
-    }
-
-    label pI = 0;
-    label seedi = events[0][0];
-
-    pointsInCell.setSize(mesh_.nCells());
-
-    const labelVector nPoints(vector::one + box_.span() / dx_);
-
-    for (label k=0; k < nPoints.z(); ++k)
-    {
-        for (label j=0; j < nPoints.y(); ++j)
-        {
-            for (label i=0; i < nPoints.x(); ++i)
-            {
-                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;
-
-                    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;
-                }
-            }
-        }
-    }
-
-    positions.shrink();
+    //- sort events by cell and in time
+    events.shrink();
 
-    weights.shrink();
+    sort(events);
 
-    for (auto& pic : pointsInCell)
-    {
-        pic.shrink();
-    }
-}
+    Info<< "Number of solidification events: "
+        << returnReduce(events.size(), sumOp<scalar>()) << endl;
 
-void Foam::foamToExaCA::write()
-{
-    if (!execute_)
-    {
-        return;
-    }
+    //- map points to cells
+    mesh_.time().cpuTimeIncrement();
+       
+    mapPoints();
+    
+    Info<< "Successfully mapped points to mesh in: "
+        << returnReduce(mesh_.time().cpuTimeIncrement(), maxOp<scalar>()) << " s"
+        << endl << endl;
 
-    // write the event data for each processor to separate files
+    //- interpolate and write ExaCA data in reduced data format
     const fileName exacaPath
     (
-        runTime_.rootPath()/runTime_.globalCaseName()/"ExaCA"
+        mesh_.time().rootPath()/mesh_.time().globalCaseName()/"ExaCA"
     );
 
     mkDir(exacaPath);
+    
+    interpolate();
+    
+    Info<< "Successfully interpolated and wrote ExaCA data in: "
+        << returnReduce(mesh_.time().cpuTimeIncrement(), maxOp<scalar>()) << " s"
+        << endl << endl;
 
-    Info<< "Number of solidification events: "
-        << returnReduce(events.size(), sumOp<scalar>()) << endl;
-
-    events.shrink();
-
-    sort(events);
-
-    // write exaca reduced data format with remelting
-    runTime_.cpuTimeIncrement();
-
-    meshSearch searchEngine(mesh_, polyMesh::CELL_TETS);
-
-    if (events.size() > 0)
-    {
-        mapPoints(searchEngine);
+    return true;
+}
 
-        interpolatePoints();
-    }
 
-    Info<< "Successfully mapped events to exaca file in: "
-        << returnReduce(runTime_.cpuTimeIncrement(), maxOp<scalar>()) << " s"
-        << endl << endl;
+bool Foam::functionObjects::ExaCA::write()
+{
+    return true;
 }
 
+
 // ************************************************************************* //
diff --git a/applications/solvers/additiveFoam/foamToExaCA/foamToExaCA.H b/applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.H
similarity index 58%
rename from applications/solvers/additiveFoam/foamToExaCA/foamToExaCA.H
rename to applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.H
index 2660446..96bfad0 100644
--- a/applications/solvers/additiveFoam/foamToExaCA/foamToExaCA.H
+++ b/applications/solvers/additiveFoam/functionObjects/ExaCA/ExaCA.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  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) 2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                 Copyright (C) 2023 Oak Ridge National Laboratory                
@@ -24,44 +24,35 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::foamToExaca
-
-Description
-    A utility for interpolating temperature data to a uniform cartestian grid
-    and writes the data in the ExaCA data format.
+    Foam::functionObjects::ExaCA
 
 SourceFiles
-    foamToExaca.C
+    ExaCA.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef foamToExaCA_H
-#define foamToExaCA_H
-    
-#include "IOdictionary.H"
-#include "volFields.H"
+#ifndef ExaCA_H
+#define ExaCA_H
+
+#include "fvMeshFunctionObject.H"
 #include "volPointInterpolation.H"
 #include "meshSearch.H"
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+namespace functionObjects
+{
 
 /*---------------------------------------------------------------------------*\
-         Class foamToExaca Declaration
+                       Class ExaCA Declaration
 \*---------------------------------------------------------------------------*/
 
-class foamToExaCA
+class ExaCA
 :
-    public IOdictionary
+    public fvMeshFunctionObject
 {
-private:
-    // Private data
-
-        const fvMesh& mesh_;
-
-        const Time& runTime_;
+    // Private Data
 
         const volScalarField& T_;
 
@@ -69,60 +60,90 @@ private:
 
         mutable pointScalarField Tp_;
 
+        DynamicList<label> overlapCells;
+
+        DynamicList<List<scalar>> events;
 
         List<DynamicList<label>> pointsInCell;
 
         DynamicList<point> positions;
+        
         DynamicList<scalarField> weights;
+              
+        boundBox  box_;
 
-        DynamicList<label> compactCells;
+        scalar    isoValue_;
 
-        DynamicList<List<scalar>> events;
+        scalar    dx_;
+        
+        meshSearch searchEngine_;
 
-        Switch    execute_;
+    // Private Member Functions
 
-        boundBox  box_;
 
-        scalar    iso_;
+public:
 
-        scalar    dx_;
+    //- Runtime type information
+    TypeName("ExaCA");
 
-public:
-    TypeName("foamToExaCA");
 
-    // Construct from temperature field
-    foamToExaCA(const volScalarField& T);
+    // Constructors
+
+        //- Construct from Time and dictionary
+        ExaCA
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+        //- Disallow default bitwise copy construction
+        ExaCA(const ExaCA&) = delete;
+
 
     //- Destructor
-    virtual ~foamToExaCA()
-    {}
+    virtual ~ExaCA();
+
 
     // Member Functions
-    virtual void initialize();
 
-    virtual void update();
+        //- Read the ExaCA data
+        virtual bool read(const dictionary&);
 
-    virtual void write();
+        //- Initialize lists
+        virtual void setOverlapCells();//initialize();
 
-    virtual void mapPoints(const meshSearch& searchEngine);
+        //- Return the list of fields required
+        virtual wordList fields() const;
 
-    virtual void interpolatePoints();
-};
+        //- Execute, currently does nothing
+        virtual bool execute();
 
+        virtual void mapPoints();
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+        virtual void interpolate();
+                                
+        //- Execute at the final time-loop, currently does nothing
+        virtual bool end();
+
+        //- Write the ExaCA data
+        virtual bool write();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const ExaCA&) = delete;
+};
 
-} // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#ifdef NoRepository
-    #include "foamToExaCA/foamToExaCA.C"
-#endif
+} // End namespace functionObjects
+} // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
 
 // ************************************************************************* //
-
diff --git a/applications/solvers/additiveFoam/functionObjects/ExaCA/Make/files b/applications/solvers/additiveFoam/functionObjects/ExaCA/Make/files
new file mode 100644
index 0000000..79b17f9
--- /dev/null
+++ b/applications/solvers/additiveFoam/functionObjects/ExaCA/Make/files
@@ -0,0 +1,3 @@
+ExaCA.C
+
+LIB = $(FOAM_USER_LIBBIN)/libExaCAFunctionObject
diff --git a/applications/solvers/additiveFoam/functionObjects/ExaCA/Make/options b/applications/solvers/additiveFoam/functionObjects/ExaCA/Make/options
new file mode 100644
index 0000000..a3ae8da
--- /dev/null
+++ b/applications/solvers/additiveFoam/functionObjects/ExaCA/Make/options
@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools
diff --git a/applications/solvers/additiveFoam/thermo/TEqn.H b/applications/solvers/additiveFoam/thermo/TEqn.H
index ad8485b..ee58df8 100644
--- a/applications/solvers/additiveFoam/thermo/TEqn.H
+++ b/applications/solvers/additiveFoam/thermo/TEqn.H
@@ -25,8 +25,12 @@
         (
             TEqn
          ==
-            rDeltaT*rho*Lf*(fvm::Sp(dFdT, T) - dFdT*T0)
-          + rho*Lf*fvc::ddt(alpha1)
+            (
+                (ddtScheme == "Euler")
+              ? rho*Lf*fv::EulerDdtScheme<scalar>(mesh).fvcDdt(alpha1)
+              : rho*Lf*fvc::ddt(alpha1)
+            )
+          + rDeltaT*rho*Lf*(fvm::Sp(dFdT, T) - dFdT*T0)
         );
 
         T.correctBoundaryConditions();
diff --git a/applications/solvers/additiveFoam/thermo/thermoScheme.H b/applications/solvers/additiveFoam/thermo/thermoScheme.H
index f004797..101a057 100644
--- a/applications/solvers/additiveFoam/thermo/thermoScheme.H
+++ b/applications/solvers/additiveFoam/thermo/thermoScheme.H
@@ -1,12 +1,30 @@
 //- Construct the temperature equation for the chosen time integration scheme
+word ddtScheme = mesh.schemes().ddt(T.name());
+
+if
+(
+    (ddtScheme != "Euler")
+ && (ddtScheme != "backward")
+ && (ddtScheme != "CrankNicolson")
+)
+{
+    FatalErrorInFunction
+        << "Unsupported ddtScheme " << ddtScheme
+        << ". Choose Euler, CrankNicolson, or backward."
+        << abort(FatalError);
+}
+
 fvScalarMatrix TEqn(T, dimPower);
 dimensionedScalar rDeltaT = 1.0 / runTime.deltaT();
-
+   
 if (explicitSolve)
 {
+    // Forward Euler is the only explicit method supported.
+    ddtScheme = "Euler";
+    
     TEqn =
     (
-        rho*Cp*(fvm::ddt(T) + fvc::div(phi, T))
+        rho*Cp*(fv::EulerDdtScheme<scalar>(mesh).fvmDdt(T) + fvc::div(phi, T))
       - fvc::laplacian(kappa, T)
       - sources.qDot()
     );
@@ -22,28 +40,15 @@ else
 
     scalar coefft = 1.0;
     
-    const word ddtScheme = mesh.schemes().ddt(T.name());
-
-    if
-    (
-        (ddtScheme != "Euler")
-     && (ddtScheme != "backward")
-     && (ddtScheme != "CrankNicolson")
-    )
-    {
-        FatalErrorInFunction
-            << "Unsupported ddtScheme " << ddtScheme
-            << ". Choose Euler, CrankNicolson, or backward."
-            << abort(FatalError);
-    }
-    else if (ddtScheme == "backward")
+    if (ddtScheme == "backward")
     {
         const dimensionedScalar deltaT = runTime.deltaT();
         const dimensionedScalar deltaT0 = runTime.deltaT0();
         
         coefft = 1.0 + (deltaT/(deltaT + deltaT0)).value();    
     }
-    else if (ddtScheme == "CrankNicolson")
+    
+    if (ddtScheme == "CrankNicolson")
     {
         scalar ocCoeff =
             refCast<const fv::CrankNicolsonDdtScheme<scalar>>
diff --git a/applications/utilities/createScanPath/Make/files b/applications/utilities/createScanPath/Make/files
new file mode 100644
index 0000000..2a874bb
--- /dev/null
+++ b/applications/utilities/createScanPath/Make/files
@@ -0,0 +1,3 @@
+createScanPath.C
+
+EXE = $(FOAM_USER_APPBIN)/createScanPath
diff --git a/applications/utilities/createScanPath/Make/options b/applications/utilities/createScanPath/Make/options
new file mode 100644
index 0000000..d27c95d
--- /dev/null
+++ b/applications/utilities/createScanPath/Make/options
@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools
diff --git a/applications/utilities/createScanPath/createFields.H b/applications/utilities/createScanPath/createFields.H
new file mode 100644
index 0000000..8e63b77
--- /dev/null
+++ b/applications/utilities/createScanPath/createFields.H
@@ -0,0 +1,384 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2024 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                Copyright (C) 2023 Oak Ridge National Laboratory                
+-------------------------------------------------------------------------------
+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/>.
+
+SourceFiles
+    createScanPath.C
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "OFstream.H"
+#include <algorithm>
+#include <fstream>
+#include <iomanip>
+#include <sstream>
+#include <vector>
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                            Struct Point Definition
+\*---------------------------------------------------------------------------*/
+
+struct Point
+{
+    scalar x;
+    scalar y;
+
+    // Default constructor
+    Point()
+        : x( 0.0 )
+        , y( 0.0 )
+    {
+    }
+
+    // Construct from x and y
+    Point( scalar xCoord, scalar yCoord )
+        : x( xCoord )
+        , y( yCoord )
+    {
+    }
+
+    // Function to rotate a point around a specified origin by a given angle
+    Point rotate( const Point& origin, scalar degrees ) const
+    {
+        scalar angle = degrees * ( M_PI / 180.0 );
+
+        scalar s = sin( angle );
+        scalar c = cos( angle );
+        scalar translatedX = x - origin.x;
+        scalar translatedY = y - origin.y;
+        scalar newX = translatedX * c - translatedY * s + origin.x;
+        scalar newY = translatedX * s + translatedY * c + origin.y;
+        return Point( newX, newY );
+    }
+
+    // Define a custom operator<< to print a Point
+    friend std::ostream& operator<<( std::ostream& os, const Point& point )
+    {
+        os << "(" << point.x << ", " << point.y << ")";
+        return os;
+    }
+};
+
+// Function to calculate the distance between two points
+scalar distance( const Point& p1, const Point& p2 )
+{
+    scalar dx = p1.x - p2.x;
+    scalar dy = p1.y - p2.y;
+    return std::sqrt( dx * dx + dy * dy );
+}
+
+
+/*---------------------------------------------------------------------------*\
+                            Struct Line Definition
+\*---------------------------------------------------------------------------*/
+
+struct Line
+{
+    Point start;
+    Point end;
+
+    // Default constructor
+    Line()
+        : start( Point() )
+        , end( Point() )
+    {
+    }
+
+    // Construct from two points
+    Line( Point startPoint, Point endPoint )
+        : start( startPoint )
+        , end( endPoint )
+    {
+    }
+
+    // Function to rotate the line around a specified origin by a given angle
+    void rotate( const Point& origin, scalar angle )
+    {
+        start = start.rotate( origin, angle );
+        end = end.rotate( origin, angle );
+    }
+
+    bool isFinite()
+    {
+        return ( !std::isnan( start.x ) && !std::isnan( start.y ) &&
+                 !std::isnan( end.x ) && !std::isnan( end.y ) );
+    }
+
+    // Define a custom operator<< to print a Line
+    friend std::ostream& operator<<( std::ostream& os, const Line& line )
+    {
+        os << "(" << line.start.x << ", " << line.start.y << "), "
+           << "(" << line.end.x << ", " << line.end.y << ")";
+        return os;
+    }
+};
+
+
+/*---------------------------------------------------------------------------*\
+                            Struct BoundBox Definition
+\*---------------------------------------------------------------------------*/
+
+struct BoundBox
+{
+    Point minPoint;
+    Point maxPoint;
+    Point midPoint;
+
+    std::vector<Line> edges;
+
+    // Construct from two points
+    BoundBox( Point minP, Point maxP )
+        : minPoint( minP )
+        , maxPoint( maxP )
+        , midPoint( ( minP.x + maxP.x ) / 2.0, ( minP.y + maxP.y ) / 2.0 )
+    {
+        // left, right, top, bottom
+        edges = { { Point( minPoint.x, minPoint.y ),
+                    Point( minPoint.x, maxPoint.y ) },
+                  { Point( maxPoint.x, minPoint.y ),
+                    Point( maxPoint.x, maxPoint.y ) },
+                  { Point( minPoint.x, maxPoint.y ),
+                    Point( maxPoint.x, maxPoint.y ) },
+                  { Point( minPoint.x, minPoint.y ),
+                    Point( maxPoint.x, minPoint.y ) } };
+    }
+
+    // Function to check if a point is inside the bounding box
+    bool isInside( Point p )
+    {
+        return p.x >= minPoint.x && p.x <= maxPoint.x && p.y >= minPoint.y &&
+               p.y <= maxPoint.y;
+    }
+
+    Line cropLine( const Line& line )
+    {
+        std::vector<Point> intersections;
+
+        for ( const Line& edge : edges )
+        {
+            Point intersection = intersect( edge, line );
+
+            if ( !std::isnan( intersection.x ) &&
+                 !std::isnan( intersection.y ) )
+            {
+                intersections.push_back( intersection );
+            }
+        }
+
+        if ( intersections.size() == 0 )
+        {
+            return Line( Point( NAN, NAN ), Point( NAN, NAN ) );
+        }
+
+        // Sort intersection points based on position along the original line
+        std::sort( intersections.begin(), intersections.end(),
+                   [&line]( const Point& p1, const Point& p2 )
+                   {
+                       scalar d1 = distance( line.start, p1 );
+                       scalar d2 = distance( line.start, p2 );
+                       return d1 < d2;
+                   } );
+
+        Line intersectedLine( intersections.front(), intersections.back() );
+
+        return intersectedLine;
+    }
+
+    // Function to find the intersection point between two lines
+    Point intersect( const Line& line1, const Line& line2 )
+    {
+        Point intersection( std::numeric_limits<scalar>::quiet_NaN(),
+                            std::numeric_limits<scalar>::quiet_NaN() );
+
+        // components of first line
+        scalar x1 = line1.start.x;
+        scalar y1 = line1.start.y;
+        scalar x2 = line1.end.x;
+        scalar y2 = line1.end.y;
+
+        // components of second line
+        scalar x3 = line2.start.x;
+        scalar y3 = line2.start.y;
+        scalar x4 = line2.end.x;
+        scalar y4 = line2.end.y;
+
+        scalar denominator =
+            ( x1 - x2 ) * ( y3 - y4 ) - ( y1 - y2 ) * ( x3 - x4 );
+
+        // Lines are parallel or colinear, no intersection
+        if ( denominator == 0 )
+        {
+            return intersection;
+        }
+
+        scalar t = ( ( x1 - x3 ) * ( y3 - y4 ) - ( y1 - y3 ) * ( x3 - x4 ) ) /
+                   denominator;
+
+        scalar u = -( ( x1 - x2 ) * ( y1 - y3 ) - ( y1 - y2 ) * ( x1 - x3 ) ) /
+                   denominator;
+
+        if ( t >= 0 && t <= 1 && u >= 0 && u <= 1 )
+        {
+            intersection.x = x1 + t * ( x2 - x1 );
+            intersection.y = y1 + t * ( y2 - y1 );
+            return intersection;
+        }
+
+        return intersection;
+    }
+};
+
+
+/*---------------------------------------------------------------------------*\
+                            Struct Path Definition
+\*---------------------------------------------------------------------------*/
+
+struct Path
+{
+    std::vector<Line> lines;
+    scalar power;
+    scalar speed;
+    scalar dwellTime;
+
+    // Construct the path from a bounding box and hatch spacing
+
+    // Default constructor
+    Path() {}
+
+    // Constructor from components
+    Path(const BoundBox& bbox, const scalar step, const scalar angle)
+    {
+        createPath(bbox, step, angle);
+    }
+    
+    // Function to create path from components
+    void createPath(BoundBox bbox, const scalar step, const scalar angle)
+    {
+        label n = countLines( bbox, step );
+
+        // create a pad of infinitely long, equally parallel lines
+        std::vector<Line> pathLines;
+
+        const scalar great = 1e10;
+
+        // Create lines in the negative direction, excluding the midpoint line
+        for ( label i = n - 1; i > 0; --i )
+        {
+            scalar height = bbox.midPoint.y - i * step;
+            Line currentLine( Point( -great, height ), Point( great, height ) );
+            pathLines.push_back( currentLine );
+        }
+
+        // Create lines in the positive direction, including the midpoint
+        for ( label i = 0; i < n; ++i )
+        {
+            scalar height = bbox.midPoint.y + i * step;
+            Line currentLine( Point( -great, height ), Point( great, height ) );
+            pathLines.push_back( currentLine );
+        }
+
+        // apply rotation and cropping to the initial lines
+        for ( const Line& pathLine : pathLines )
+        {
+            // Rotate the endpoints by the specified angle
+            Line rotatedLine = pathLine;
+            rotatedLine.rotate( bbox.midPoint, angle );
+
+            Line croppedLine = bbox.cropLine( rotatedLine );
+
+            if ( croppedLine.isFinite() )
+            {
+                lines.push_back( croppedLine );
+            }
+        }
+    }
+
+    // Function to find the number of scan vectors in the bounding box
+    label countLines( const BoundBox& bbox, scalar step ) const
+    {
+        label nX = 0;
+        label nY = 0;
+
+        for ( scalar x = bbox.minPoint.x; x <= bbox.maxPoint.x; x += step )
+        {
+            nX++;
+        }
+
+        for ( scalar y = bbox.minPoint.y; y <= bbox.maxPoint.y; y += step )
+        {
+            nY++;
+        }
+
+        return ( nX > nY ) ? nX : nY;
+    }
+
+    void write( const std::string& filename,
+                const bool bi_direction = true ) const
+    {
+        std::ofstream file( filename );
+
+        file << "Mode\tX(m)\tY(m)\tZ(m)\tPower(W)\ttParam" << std::endl;
+
+        for ( size_t i = 0; i < lines.size(); ++i )
+        {
+            const Line& line = lines[i];
+
+            Point first = line.start;
+            Point second = line.end;
+
+            // reverse the odd lines for bi_directional
+            if ( bi_direction && i % 2 == 1 )
+            {
+                first = line.end;
+                second = line.start;
+            }
+
+            // hatch (with skywrite)
+            if ( i == 0 )
+            {
+                // no initial dwell
+                file << "1\t" << first.x << "\t" << first.y << "\t0\t0\t0\n";
+            }
+            else
+            {
+                file << "1\t" << first.x << "\t" << first.y << "\t0\t" << 0
+                     << "\t" << dwellTime << "\n";
+            }
+
+            // raster
+            file << "0\t" << second.x << "\t" << second.y << "\t0\t" << power
+                 << "\t" << speed << "\n";
+        }
+
+        file.close();
+    }
+};
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/applications/utilities/createScanPath/createScanPath.C b/applications/utilities/createScanPath/createScanPath.C
new file mode 100644
index 0000000..9667dfc
--- /dev/null
+++ b/applications/utilities/createScanPath/createScanPath.C
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2024 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+                Copyright (C) 2023 Oak Ridge National Laboratory                
+-------------------------------------------------------------------------------
+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
+    createScanPath
+
+\*---------------------------------------------------------------------------*/
+
+#include "createFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    
+    IOdictionary dict
+    (
+        IOobject
+        (
+            "createScanPathDict",
+            runTime.constant(),
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+    
+    point2D p0 = dict.lookup<point2D>("minPoint");
+    Point minPoint(p0.x(), p0.y());
+
+    point2D p1 = dict.lookup<point2D>("maxPoint");
+    Point maxPoint(p1.x(), p1.y());   
+    
+    scalar angle = dict.lookup<scalar>("angle");
+    scalar hatch = dict.lookup<scalar>("hatch");
+    
+    label nRotations = dict.lookup<label>("nRotations");
+    
+    scalar power = dict.lookup<scalar>("power");
+    scalar speed = dict.lookup<scalar>("speed");
+    scalar dwellTime = dict.lookup<scalar>("dwellTime");
+    
+    bool biDirection = dict.lookupOrDefault<bool>("biDirection", true);
+    
+    // Create bounding box for scan vectors
+    BoundBox boundingBox( minPoint, maxPoint );
+
+    // Rotate and write scan vectors to file
+    scalar rotation = 0.0;
+    
+    for (label i = 0; i < nRotations; ++i)
+    {
+        Info << "Creating scan path for rotation angle: " << rotation << endl;
+        
+        Path path( boundingBox, hatch, rotation );
+        
+        path.power = power;
+        path.speed = speed;
+        path.dwellTime = dwellTime;
+
+        std::ostringstream oss;
+        oss << std::fixed << std::setprecision(0) << rotation;
+        std::string rotationString = oss.str();
+        
+        const fileName filename
+        (
+            runTime.constant()/"scanPath_" + std::to_string(i)
+        );
+
+        path.write(filename, biDirection);
+
+        rotation += angle;
+    }
+
+    Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/AMB2018-02-B/constant/foamToExaCADict b/applications/utilities/createScanPath/createScanPathDict
similarity index 75%
rename from tutorials/AMB2018-02-B/constant/foamToExaCADict
rename to applications/utilities/createScanPath/createScanPathDict
index c874ff4..405dfee 100644
--- a/tutorials/AMB2018-02-B/constant/foamToExaCADict
+++ b/applications/utilities/createScanPath/createScanPathDict
@@ -10,15 +10,27 @@ FoamFile
     version     2.0;
     format      ascii;
     class       dictionary;
-    object      foamToExaCADict;
+    object      createScanPathDict;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-execute     on;
+minPoint    (0.0 0.0);
 
-box         (0 -0.0001 -0.0002) (0.002 0.0001 0);
-dx          5e-6;
-isotherm    1620;
+maxPoint    (0.002 0.002);
+
+angle       67;
+
+hatch       1e-4;
+
+nRotations  5;
+
+power       200;
+
+speed       1;
+
+dwellTime   0.0;
+
+biDirection true;
 
 // ************************************************************************* //
diff --git a/applications/utilities/multiLayer/runLayers b/applications/utilities/multiLayer/runLayers
index 7877ce5..17552dd 100755
--- a/applications/utilities/multiLayer/runLayers
+++ b/applications/utilities/multiLayer/runLayers
@@ -153,9 +153,23 @@ do
     fi
 
     thickness="$(foamDictionary -entry linearDirectionCoeffs/thickness -value system/extrudeMeshDict)"
+    
+    # set powder height to the layer thickness 
+    foamDictionary -entry regions -set \
+    "(
+        boxToCell
+        {
+            box (-1 -1 -$layerThickness) (1 1 1);
+            fieldValues
+            (
+                volScalarFieldValue alpha.powder 1
+            );
+        }
+    )" system/setFieldsDict
+    
     nProcs="$(foamDictionary -entry numberOfSubdomains -value system/decomposeParDict)"
 
-    transformPoints "translate=(0 0 -$thickness)"
+    transformPoints "translate=(0 0 -$thickness)"    
     setFields
     decomposePar
     mpirun -np $nProcs additiveFoam -parallel
diff --git a/doc/paper.bib b/doc/paper.bib
new file mode 100644
index 0000000..32e0070
--- /dev/null
+++ b/doc/paper.bib
@@ -0,0 +1,193 @@
+@article{openfoam,
+   abstract = {In this article the principles of the field operation and manipulation (FOAM) C++class library for continuum mechanics are outlined. Our intention is to make it as easy as possible to develop reliable and efficient computational continuum-mechanics codes: this is achieved by making the top-level syntax of the code as close as possible to conventional mathematical notation for tensors and partial differential equations. Object-orientation techniques enable the creation of data types that closely mimic those of continuum mechanics, and the operator overloading possible in C++ allows normal mathematical symbols to be used for the basic operations. As an example, the implementation of various types of turbulence modeling in a FOAM computational-fluid-dynamics code is discussed, and calculations performed on a standard test case, that of flow around a square prism, are presented. To demonstrate the flexibility of the FOAM library, codes for solving structures and magnetohydrodynamics are also presented with appropriate test case results given.},
+   author = {H. G. Weller and G. Tabor and H. Jasak and C. Fureby},
+   doi = {10.1063/1.168744},
+   issn = {08941866},
+   issue = {6},
+   journal = {Computers in Physics},
+   pages = {620},
+   title = {A tensorial approach to computational continuum mechanics using object-oriented techniques},
+   volume = {12},
+   year = {1998},
+}
+
+@techreport{ALE3D,
+  title={ALE3D: An arbitrary Lagrangian-Eulerian multi-physics code},
+  author={Noble, Charles R and Anderson, Andrew T and Barton, Nathan R and Bramwell, Jamie A and Capps, Arlie and Chang, Michael H and Chou, Jin J and Dawson, David M and Diana, Emily R and Dunn, Timothy A and others},
+  year={2017},
+  institution={Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States)}
+}
+
+@Manual{FLOW-3D,
+  author = {Flow Science, Inc.},
+  title = {FLOW-3D, Version~2023R1},
+  year = {2023},
+  address = {Santa Fe, NM},
+  url = {https://www.flow3d.com/}
+}
+
+@manual{abaqus,
+  title        = {Abaqus Documentation},
+  author       = {{Dassault Systèmes Simulia Corp.}},
+  organization = {Dassault Systèmes},
+  address      = {Providence, RI, USA},
+  year         = {2024},
+  note         = {Version 2024},
+  url          = {https://www.3ds.com/products-services/simulia/products/abaqus/}
+}
+
+@manual{ansys,
+  title        = {ANSYS Documentation},
+  author       = {{ANSYS Inc.}},
+  organization = {ANSYS Inc.},
+  address      = {Canonsburg, PA, USA},
+  year         = {2024},
+  note         = {Release 2024 R1},
+  url          = {https://www.ansys.com/products/structures/ansys-mechanical}
+}
+
+@software{adamantine,
+  author = {Turcksin, Bruno and DeWitt, Stephen},
+  doi = {10.11578/dc.20211129.1},
+  title = {{adamantine}},
+  url = {https://github.com/adamantine-sim/adamantine}
+}
+
+@Article{deal.ii,
+  title   = {The \texttt{deal.II} Library, Version 9.5},
+  author  = {Daniel Arndt and Wolfgang Bangerth and Maximilian Bergbauer and
+             Marco Feder and Marc Fehling and Johannes Heinz and
+             Timo Heister and Luca Heltai and Martin Kronbichler and
+             Matthias Maier and Peter Munch and Jean-Paul Pelteret and
+             Bruno Turcksin and David Wells and Stefano Zampini},
+  journal = {Journal of Numerical Mathematics},
+  year    = {2023},
+  doi     = {10.1515/jnma-2023-0089},
+  pages   = {231--246},
+  volume  = {31},
+  number  = {3},
+  url     = {https://dealii.org/deal95-preprint.pdf}
+}
+
+@article{exaca,
+  title={ExaCA: A performance portable exascale cellular automata application for alloy solidification modeling},
+  author={Rolchigo, Matt and Reeve, Samuel Temple and Stump, Benjamin and Knapp, Gerald L and Coleman, John and Plotkowski, Alex and Belak, James},
+  journal={Computational Materials Science},
+  volume={214},
+  pages={111692},
+  year={2022},
+  publisher={Elsevier}
+}
+
+@article{texture-regimes,
+title = {Grain structure and texture selection regimes in metal powder bed fusion},
+journal = {Additive Manufacturing},
+volume = {81},
+pages = {104024},
+year = {2024},
+issn = {2214-8604},
+doi = {10.1016/j.addma.2024.104024},
+url = {https://www.sciencedirect.com/science/article/pii/S2214860424000708},
+author = {Matt Rolchigo and John Coleman and Gerry L. Knapp and Alex Plotkowski}
+}
+
+@article{knapp,
+title = {Calibrating uncertain parameters in melt pool simulations of additive manufacturing},
+journal = {Computational Materials Science},
+volume = {218},
+pages = {111904},
+year = {2023},
+issn = {0927-0256},
+doi = {10.1016/j.commatsci.2022.111904},
+url = {https://www.sciencedirect.com/science/article/pii/S0927025622006152},
+author = {G.L. Knapp and J. Coleman and M. Rolchigo and M. Stoyanov and A. Plotkowski}
+}
+
+@article{haines-1,
+  title = {Experimental and computational analysis of site-specific formation of phases in laser powder bed fusion 17–4 precipitate hardened stainless steel},
+  journal = {Additive Manufacturing},
+  volume = {73},
+  pages = {103686},
+  year = {2023},
+  issn = {2214-8604},
+  doi = {10.1016/j.addma.2023.103686},
+  url = {https://www.sciencedirect.com/science/article/pii/S2214860423002993},
+  author = {Michael P. Haines and Maxwell S. Moyle and Vitor V. Rielli and Vladimir Luzin and Nima Haghdadi and Sophie Primig}
+}
+
+@article{haines-2,
+  title = {Role of scan strategies and heat treatment on grain structure evolution in Fe-Si soft magnetic alloys made by laser-powder bed fusion},
+  journal = {Additive Manufacturing},
+  volume = {50},
+  pages = {102578},
+  year = {2022},
+  issn = {2214-8604},
+  doi = {10.1016/j.addma.2021.102578},
+  url = {https://www.sciencedirect.com/science/article/pii/S2214860421007259},
+  author = {MP Haines and F. List and K. Carver and DN Leonard and A. Plotkowski and CM Fancher and RR Dehoff and SS Babu}
+}
+
+@article{plotkowski,
+  author = {Plotkowski, A. and Coleman, J. and Fancher, C. M. and Haines, M. P. and Whetten, S. R. and Kustas, A. B.},
+  year = {2024},
+  title = {Grain Structure Evolution in Fe-6Si During Directed Energy Deposition},
+  journal = {JOM},
+  pages = {1031},
+  volume = {76},
+  number = {3},
+  doi = {10.1007/s11837-023-06279-3}
+}
+
+@software{AdditiveFOAM_1.0.0,
+  author       = {John Coleman and
+                  Kellis Kincaid and
+                  Gerald L. Knapp and
+                  Benjamin Stump and
+                  Alexander J. Plotkowski},
+  title        = {AdditiveFOAM: Release 1.0},
+  month        = jun,
+  year         = 2023,
+  publisher    = {Zenodo},
+  version      = {1.0.0},
+  doi          = {10.5281/zenodo.8034098}
+}
+
+@INPROCEEDINGS{zoltan,
+ AUTHOR = "K.D. Devine and E.G. Boman and L.A. Riesen and U.V. Catalyurek 
+          and C. Chevalier",
+ TITLE = "Getting Started with Zoltan:  A Short Tutorial",
+ BOOKTITLE = "Proc. of 2009 Dagstuhl Seminar on Combinatorial Scientific 
+              Computing", 
+ YEAR = "2009",
+ NOTE = "Also available as Sandia National Labs Tech Report SAND2009-0578C."
+}
+
+@article{amb2018,
+  title={Outcomes and conclusions from the 2018 AM-Bench measurements, challenge problems, modeling submissions, and conference},
+  author={Levine, Lyle E. and
+          Lane, Brandon M. and
+          Heigel, Jarred C. and
+          Migler, Kalman and
+          Stoudt, Mark R. and
+          Phan, Thien Q. and
+          Ricker, Richard E. and
+          Strantza, Maria and
+          Hill, Martin R. and
+          Zhang, Fan and
+          Seppala, Jonathan and
+          Garboczi, Edward and
+          Bain, Erich and
+          Cole, Daniel and
+          Allen, Andrew and
+          Fox, Jason and
+          Campbell, Carelyn},
+  journal={Integrating Materials and Manufacturing Innovation},
+  volume={9},
+  number={3},
+  pages={362--377},
+  year={2020},
+  publisher={Springer},
+  doi={10.1007/s40192-019-00164-1}
+}
+
+
diff --git a/doc/paper.md b/doc/paper.md
new file mode 100644
index 0000000..5fe8d83
--- /dev/null
+++ b/doc/paper.md
@@ -0,0 +1,75 @@
+---
+title: 'AdditiveFOAM: A Continuum Multiphysics Code for Additive Manufacturing'
+tags:
+  - Computational Fluid Dynamics
+  - Heat Transfer
+  - Mass Transfer
+  - Additive Manufacturing
+  - OpenFOAM
+authors:
+  - name: John Coleman
+    affiliation: "1"
+    corresponding: true
+  - name: Kellis Kincaid
+    affiliation: "2"
+  - name: Gerry L. Knapp
+    affiliation: "3"
+  - name: Benjamin Stump
+    affiliation: "1"
+  - name: Samuel Reeve
+    affiliation: "1"
+  - name: Matt Rolchigo
+    affiliation: "1"
+  - name: Alex Plotkowski
+    affiliation: "1"
+affiliations:
+ - name: Computational Science and Engineering Division, Oak Ridge National Laboratory
+   index: 1
+ - name: Nuclear Energy and Fuel Cycle Division, Oak Ridge National Laboratory
+   index: 2
+ - name: Material Science and Technology Division, Oak Ridge National Laboratory
+   index: 3
+date: 12 June 2024
+bibliography: paper.bib
+---
+
+# Summary
+
+AdditiveFOAM is a computational framework that simulates transport phenomena in Additive Manufacturing (AM) processes. It is built on OpenFOAM [@openfoam], the leading free, open-source software package for computational fluid dynamics (CFD). OpenFOAM offers an extensible platform for solving complex multiphysics problems using state-of-the-art finite volume methods. AdditiveFOAM leverages these capabilities to develop specialized tools aimed at addressing challenges in AM processing.
+
+Metal additive manufacturing, also known as metal 3D printing, is an advanced manufacturing technique that creates physical parts from a three-dimensional (3D) digital model by melting metal powder or wire feedstock. A significant area of research in metal AM focuses on process planning to mitigate anomalous features during printing that are deleterious to part performance (e.g., porosity and cracking), as well as controlling localized microstructure and material properties. Given the high costs and substantial time requirements associated with experimental methods for qualifying new materials and processes, there is a compelling incentive for researchers to utilize advanced computational simulations. In this context, AdditiveFOAM offers a simulation framework to better understand undesirable features in printing, thereby enhancing process planning and reducing the reliance on labor-intensive experimental campaigns.
+
+# Statement of need
+
+Numerical models for AM traditionally seek to represent physics at a single length scale depending on the target question. Two main categories of models with available software solutions are distinguished: 1) high-fidelity melt pool-scale models, and 2) part-scale heat transfer models.
+
+The first category contains models that seek to resolve all relevant physics within the melt pool. These models can produce excellent agreement with experiments and offer insights into the underlying causes of anomalous features during processing, but are generally expensive to run, limiting the scan length that can be simulated to a few millimeters with current computational resources. Due to their computational expense, these models are well-suited for use at specialized research institutions with large HPC infrastructure to answer target questions about melt pool scale physical phenomena. For example, ALE3D [@ALE3D] is a versatile multiphysics simulation tool that uses the Arbitrary Lagrangian-Eulerian approach and has been used for powder-resolved simulations of laser-material interactions in AM. However, ALE3D is a limited access code for use by United States Department of Defense/Energy laboratories and their contractors, while ALE3D4I is available for U.S. companies and academics through individual use agreements. Alternatively, FLOW-3D [@FLOW-3D] is a commercial CFD software known for its capabilities in simulating complex free-surface problems, and it has several specialized models for AM. However, FLOW-3D is proprietary, meaning its source code is not available for public inspection or modification, and users must purchase a license to use the software.
+
+The second category contains models that simplify melt pool physics to rapidly simulate heat transport, residual stress, and distortion across an entire part. Commercial finite element thermomechanics software solutions built on Abaqus [@abaqus] and Ansys [@ansys] are available; however, these software packages are not open and free to use and develop upon. Alternatively, Adamantine [@adamantine] is a thermomechanics simulation tool built on an open-source software stack designed for high-performance computing across various architectures, including the deal.II finite element software package [@deal.ii], providing an alternative to proprietary software solutions.
+
+There is an established need for intermediate frameworks that accurately predict melt pool shape, thermal gradients, and other meso-scale quantities across an entire part. These models can inform process design decisions through heuristic estimations of anomalous printing features (e.g., keyhole formation and lack-of-fusion) and predict the final microstructure and material properties of AM components. AdditiveFOAM addresses these challenges, featuring a volumetric source term in the energy equation for multiple heat sources, optional Marangoni-driven fluid flow, and a scheme for implementing tabulated thermodynamic pathways for metal alloys.
+
+# Software Features
+The main feature of AdditiveFOAM is a transient, multiphysics application that solves conservation equations for mass, momentum, and energy during AM processing, built upon the OpenFOAM finite volume software package [@openfoam]. This solver includes a novel thermodynamic algorithm enabling variable-order time integration for tabulated thermodynamic pathways. The available time integration schemes are explicit forward Euler, implicit backward Euler, backward differentiation formula (BDF-2), and Crank-Nicolson. An adaptive time integration and time-stepping method automatically switches between implicit and explicit schemes to balance computational cost and solution accuracy depending on the problem state and user-defined tolerances. Additionally, AdditiveFOAM features boundary conditions for Marangoni-driven fluid flow in the melt pool, as well as convective and radiative heat transfer. Another feature of AdditiveFOAM is the volumetric heat source library that supports any number of independently moving sources with distinct energy profiles that have been validated for several AM processes, including laser powder bed fusion (LPBF), directed energy deposition (DED), and electron beam melting (EBM). Finally, AdditiveFOAM leverages existing OpenFOAM libraries to sample thermal data needed for microstructure predictions with a specific library for coupling to the grain structure prediction software ExaCA [@exaca]. Future releases of AdditiveFOAM will focus on resource-optimized adaptive mesh refinement (AMR), dynamic load-balancing using the Zoltan [@zoltan] library, and the integration of next-generation laser profiles (e.g., nLight AFX series) into the volumetric heat source library.
+
+AdditiveFOAM was designed to be used openly by researchers, industry, and academics. It has already been used in several scientific publications primarily focused on studying how AM process planning influences the development of local microstructural features in a part. Select publications that demonstrate AdditiveFOAM’s usage towards understanding the connection between process planning and microstructure evolution in metal AM include:
+
+- Knapp et al. [@knapp] developed an automated framework for statistical calibration of the gaussian heat source parameters available in AdditiveFOAM to accurately predict the melt pool shape, solidification conditions, and solidification grain structure for single track welds on IN625 plate. 
+
+- Rolchigo et al. [@exaca; @texture-regimes] used AdditiveFOAM to generate thermal data for solidification grain structure predictions is laser powder bed fusion, validated against the NIST 2018 AM-Bench dataset [@amb2018]. 
+
+- Haines et al. used AdditiveFOAM to confirm experimentally observed phase transformation pathways in 17–4 PH stainless steel during the laser powder bed fusion (L-PBF) process [@haines-1], and investigated the role of scan strategy on recrystallization and grain growth in Fe-Si steel made by L-PBF with a pulsed laser [@haines-2].
+
+- Plotkowski et al. [@plotkowski] used AdditiveFOAM to confirm experimentally observed grain growth in Fe-Si steel during the Laser Engineered Net Shaping (LENS) process. 
+
+The free, open-source nature of AdditiveFOAM will continue to enable scientific explorations in AM processing. Documentation for AdditiveFOAM is hosted on GitHub pages with the link: https://ornl.github.io/AdditiveFOAM/. 
+
+# Acknowledgements
+
+This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy (DOE). The publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan. 
+
+This research used resources of the Oak Ridge Leadership Computing Facility (OLCF), supported by DOE under contract DE-AC05-00OR22725. 
+
+This research used resources of the Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. 
+
+# References
diff --git a/tutorials/AMB2018-02-B/Allclean b/tutorials/AMB2018-02-B/Allclean
index 6b45f82..5c72b39 100755
--- a/tutorials/AMB2018-02-B/Allclean
+++ b/tutorials/AMB2018-02-B/Allclean
@@ -5,5 +5,7 @@ cd ${0%/*} || exit 1    # Run from this directory
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 cleanCase
-
+#------------------------------------------------------------------------------
+rm -rf ExaCA/Output*
+rm -rf ExaCA/time-temperature.csv
 #------------------------------------------------------------------------------
diff --git a/tutorials/AMB2018-02-B/Allrun b/tutorials/AMB2018-02-B/Allrun
index 9c9d763..ba6eb2c 100755
--- a/tutorials/AMB2018-02-B/Allrun
+++ b/tutorials/AMB2018-02-B/Allrun
@@ -3,9 +3,27 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
-
 application=`getApplication`
 
+# Parse arguments
+withExaCA=false
+while [ "$#" -gt 0 ]; do
+  case "$1" in
+    -withExaCA)
+      withExaCA=true
+      ;;
+  esac
+  shift
+done
+
+if [ "$withExaCA" = true ]; then
+    foamDictionary -entry functions -set "{#includeFunc ExaCA}" system/controlDict
+else
+    foamDictionary -entry functions -set "{}" system/controlDict
+fi
+
+#------------------------------------------------------------------------------
+# AdditiveFOAM
 runApplication blockMesh
 runApplication decomposePar
 runParallel $application
@@ -13,3 +31,13 @@ runParallel $application
 runApplication reconstructPar
 
 #------------------------------------------------------------------------------
+# ExaCA (Version 2.0.1)
+if [ "$withExaCA" = true ]; then
+    echo "x,y,z,tm,ts,cr" > ExaCA/time-temperature.csv
+    cat ExaCA/data_* >> ExaCA/time-temperature.csv
+    rm -rf ExaCA/data_*
+
+    NPROCS="$(foamDictionary -expand -entry numberOfSubdomains -value system/decomposeParDict)"
+    mpirun -np $NPROCS ~/install/exaca/bin/ExaCA ExaCA/input.json
+fi
+#------------------------------------------------------------------------------
diff --git a/tutorials/AMB2018-02-B/ExaCA/input.json b/tutorials/AMB2018-02-B/ExaCA/input.json
new file mode 100644
index 0000000..f0676a2
--- /dev/null
+++ b/tutorials/AMB2018-02-B/ExaCA/input.json
@@ -0,0 +1,35 @@
+{
+   "SimulationType": "FromFile",
+   "MaterialFileName": "Inconel625.json",
+   "GrainOrientationFile": "GrainOrientationVectors.csv",
+   "RandomSeed": 0,
+   "Domain": {
+       "CellSize": 2.5,
+       "TimeStep": 0.125,
+       "NumberOfLayers": 1,
+       "LayerOffset": 0
+   },
+   "Nucleation": {
+      "Density": 100,
+      "MeanUndercooling": 5,
+      "StDev": 0.5
+   },
+   "TemperatureData": {
+       "LayerwiseTempRead": false,
+       "TemperatureFiles": ["./ExaCA/time-temperature.csv"]
+   },
+   "Substrate": {
+      "MeanBaseplateGrainSize": 8.3,
+      "MeanPowderGrainSize": 2.5
+   },
+   "Printing": {
+      "PathToOutput": "./ExaCA/",
+      "OutputFile": "Output",
+      "PrintBinary": false,
+      "PrintExaConstitSize": 0,
+      "Interlayer": {
+          "Fields": ["GrainID", "LayerID", "GrainMisorientation", "UndercoolingCurrent", "SolidificationEventCounter"],
+          "Layers": [0, 1]
+      }
+   }
+}
diff --git a/tutorials/AMB2018-02-B/constant/heatSourceDict b/tutorials/AMB2018-02-B/constant/heatSourceDict
index 0122d92..82342cb 100644
--- a/tutorials/AMB2018-02-B/constant/heatSourceDict
+++ b/tutorials/AMB2018-02-B/constant/heatSourceDict
@@ -17,6 +17,12 @@ FoamFile
 
 sources (beam);
 
+/*---------------------------------------------------------------------------*\
+Citation:
+    G.L. Knapp, J. Coleman, M. Rolchigo, M. Stoyanov, A. Plotkowski,
+    Calibrating uncertain parameters in melt pool simulations of additive 
+    manufacturing (2023), https://doi.org/10.1016/j.commatsci.2022.11190.
+\*---------------------------------------------------------------------------*/
 beam
 {
     pathName            scanPath;    
@@ -25,7 +31,7 @@ beam
     
     constantCoeffs
     {
-        eta             0.3;
+        eta             0.33;
     }
     
     heatSourceModel     superGaussian;
diff --git a/tutorials/marangoni/constant/foamToExaCADict b/tutorials/AMB2018-02-B/system/ExaCA
similarity index 59%
rename from tutorials/marangoni/constant/foamToExaCADict
rename to tutorials/AMB2018-02-B/system/ExaCA
index 0498e97..fb9b4c4 100644
--- a/tutorials/marangoni/constant/foamToExaCADict
+++ b/tutorials/AMB2018-02-B/system/ExaCA
@@ -5,20 +5,15 @@
     \\  /    A nd           | Version:  10
      \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
-FoamFile
+
+ExaCA
 {
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      foamToExaCADict;
+    type ExaCA;
+    libs ("libExaCAFunctionObject.so");
+    
+    box         (0 -0.0001 -0.0002) (0.002 0.0001 0);
+    dx          2.5e-6;
+    isoValue    1620;
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-execute     off;
-
-box         (0 -0.0001 -0.0002) (0.002 0.0001 0);
-dx          5e-6;
-isotherm    1620;
-
-// ************************************************************************* //
+// *********************************************************************** //
diff --git a/tutorials/AMB2018-02-B/system/controlDict b/tutorials/AMB2018-02-B/system/controlDict
index ad97bac..6045cde 100644
--- a/tutorials/AMB2018-02-B/system/controlDict
+++ b/tutorials/AMB2018-02-B/system/controlDict
@@ -7,11 +7,11 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      controlDict;
+    version         2;
+    format          ascii;
+    class           dictionary;
+    location        "system";
+    object          controlDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -25,7 +25,7 @@ stopAt          endTime;
 
 endTime         0.004;
 
-deltaT          1e-7;
+deltaT          1e-07;
 
 writeControl    adjustableRunTime;
 
@@ -48,6 +48,15 @@ runTimeModifiable yes;
 adjustTimeStep  yes;
 
 maxCo           0.5;
+
 maxDi           1;
+
 maxAlphaCo      1;
+
+functions
+{
+    #includeFunc    ExaCA
+}
+
+
 // ************************************************************************* //
diff --git a/tutorials/multiBeam/Allclean b/tutorials/multiBeam/Allclean
index 6b45f82..5c72b39 100755
--- a/tutorials/multiBeam/Allclean
+++ b/tutorials/multiBeam/Allclean
@@ -5,5 +5,7 @@ cd ${0%/*} || exit 1    # Run from this directory
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 cleanCase
-
+#------------------------------------------------------------------------------
+rm -rf ExaCA/Output*
+rm -rf ExaCA/time-temperature.csv
 #------------------------------------------------------------------------------
diff --git a/tutorials/multiBeam/Allrun b/tutorials/multiBeam/Allrun
index 9c9d763..ba6eb2c 100755
--- a/tutorials/multiBeam/Allrun
+++ b/tutorials/multiBeam/Allrun
@@ -3,9 +3,27 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 # Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
-
 application=`getApplication`
 
+# Parse arguments
+withExaCA=false
+while [ "$#" -gt 0 ]; do
+  case "$1" in
+    -withExaCA)
+      withExaCA=true
+      ;;
+  esac
+  shift
+done
+
+if [ "$withExaCA" = true ]; then
+    foamDictionary -entry functions -set "{#includeFunc ExaCA}" system/controlDict
+else
+    foamDictionary -entry functions -set "{}" system/controlDict
+fi
+
+#------------------------------------------------------------------------------
+# AdditiveFOAM
 runApplication blockMesh
 runApplication decomposePar
 runParallel $application
@@ -13,3 +31,13 @@ runParallel $application
 runApplication reconstructPar
 
 #------------------------------------------------------------------------------
+# ExaCA (Version 2.0.1)
+if [ "$withExaCA" = true ]; then
+    echo "x,y,z,tm,ts,cr" > ExaCA/time-temperature.csv
+    cat ExaCA/data_* >> ExaCA/time-temperature.csv
+    rm -rf ExaCA/data_*
+
+    NPROCS="$(foamDictionary -expand -entry numberOfSubdomains -value system/decomposeParDict)"
+    mpirun -np $NPROCS ~/install/exaca/bin/ExaCA ExaCA/input.json
+fi
+#------------------------------------------------------------------------------
diff --git a/tutorials/multiBeam/ExaCA/input.json b/tutorials/multiBeam/ExaCA/input.json
new file mode 100644
index 0000000..f0676a2
--- /dev/null
+++ b/tutorials/multiBeam/ExaCA/input.json
@@ -0,0 +1,35 @@
+{
+   "SimulationType": "FromFile",
+   "MaterialFileName": "Inconel625.json",
+   "GrainOrientationFile": "GrainOrientationVectors.csv",
+   "RandomSeed": 0,
+   "Domain": {
+       "CellSize": 2.5,
+       "TimeStep": 0.125,
+       "NumberOfLayers": 1,
+       "LayerOffset": 0
+   },
+   "Nucleation": {
+      "Density": 100,
+      "MeanUndercooling": 5,
+      "StDev": 0.5
+   },
+   "TemperatureData": {
+       "LayerwiseTempRead": false,
+       "TemperatureFiles": ["./ExaCA/time-temperature.csv"]
+   },
+   "Substrate": {
+      "MeanBaseplateGrainSize": 8.3,
+      "MeanPowderGrainSize": 2.5
+   },
+   "Printing": {
+      "PathToOutput": "./ExaCA/",
+      "OutputFile": "Output",
+      "PrintBinary": false,
+      "PrintExaConstitSize": 0,
+      "Interlayer": {
+          "Fields": ["GrainID", "LayerID", "GrainMisorientation", "UndercoolingCurrent", "SolidificationEventCounter"],
+          "Layers": [0, 1]
+      }
+   }
+}
diff --git a/tutorials/multiBeam/constant/heatSourceDict b/tutorials/multiBeam/constant/heatSourceDict
index a61749a..6c7690b 100644
--- a/tutorials/multiBeam/constant/heatSourceDict
+++ b/tutorials/multiBeam/constant/heatSourceDict
@@ -17,7 +17,12 @@ FoamFile
 
 sources (beam1 beam2);
 
-//- static heat source with a gaussian distribution
+/*---------------------------------------------------------------------------*\
+Citation:
+    G.L. Knapp, J. Coleman, M. Rolchigo, M. Stoyanov, A. Plotkowski,
+    Calibrating uncertain parameters in melt pool simulations of additive 
+    manufacturing (2023), https://doi.org/10.1016/j.commatsci.2022.11190.
+\*---------------------------------------------------------------------------*/
 beam1
 {
     pathName            scanPath_1;
@@ -25,7 +30,7 @@ beam1
     absorptionModel     constant;    
     constantCoeffs
     {
-        eta             0.35;
+        eta             0.33;
     }   
 
     heatSourceModel     superGaussian;
@@ -38,7 +43,12 @@ beam1
     }
 }
 
-//- transient heat source with a modified super-gaussian distribution
+/*---------------------------------------------------------------------------*\
+Citation:
+    J. Coleman, G.L. Knapp, B. Stump, M. Rolchigo, K. Kincaid, A. Plotkowski,
+    A dynamic volumetric heat source model for laser additive manufacturing,
+    Additive Manufacturing (2024), https://doi.org/10.1016/j.addma.2024.104531.
+\*---------------------------------------------------------------------------*/
 beam2
 {
     pathName            scanPath_2;
@@ -46,19 +56,19 @@ beam2
     absorptionModel     Kelly;    
     KellyCoeffs
     {
-        geometry        cylinder;
-        eta0            0.3;
+        geometry        cone;
+        eta0            0.28;
         etaMin          0.35;
     }   
 
     heatSourceModel     modifiedSuperGaussian;
     modifiedSuperGaussianCoeffs
     {
-        dimensions      (50.0e-6 50.0e-6 30e-6);
-        m               2.0;
-        k               2.0;
+        dimensions      (40.0e-6 40.0e-6 30e-6);
+        m               2.72;
+        k               7.95;
         transient       true;
-        isoValue        1410;
+        isoValue        1620;
         nPoints         (10 10 10);
     }
 }
diff --git a/tutorials/multiBeam/constant/scanPath_1 b/tutorials/multiBeam/constant/scanPath_1
index 860d64e..9b92b66 100644
--- a/tutorials/multiBeam/constant/scanPath_1
+++ b/tutorials/multiBeam/constant/scanPath_1
@@ -1,3 +1,3 @@
 Mode    X       Y       Z   Power   Param
-1       0.000   -10e-5   0   0       0
-0       0.002   -10e-5   0   179.2   0.8
+1       0.000   -5e-5   0   0       0
+0       0.002   -5e-5   0   179.2   0.8
diff --git a/tutorials/multiBeam/constant/scanPath_2 b/tutorials/multiBeam/constant/scanPath_2
index 329fedc..9823096 100644
--- a/tutorials/multiBeam/constant/scanPath_2
+++ b/tutorials/multiBeam/constant/scanPath_2
@@ -1,3 +1,3 @@
 Mode    X       Y       Z   Power   Param
-1       0.002   10e-5   0   0       0
-0       0.000   10e-5   0   179.2   0.8
+1       0.002   5e-5   0   0       0
+0       0.000   5e-5   0   195.0   0.8
diff --git a/tutorials/multiBeam/system/ExaCA b/tutorials/multiBeam/system/ExaCA
new file mode 100644
index 0000000..83561cd
--- /dev/null
+++ b/tutorials/multiBeam/system/ExaCA
@@ -0,0 +1,19 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  10
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+
+ExaCA
+{
+    type ExaCA;
+    libs ("libExaCAFunctionObject.so");
+    
+    box         (0.0005 -0.00025 -0.0002) (0.0015 0.00025 0);
+    dx          2.5e-6;
+    isoValue    1620;
+}
+
+// *********************************************************************** //
diff --git a/tutorials/multiBeam/system/blockMeshDict b/tutorials/multiBeam/system/blockMeshDict
index c96b440..854f5ad 100644
--- a/tutorials/multiBeam/system/blockMeshDict
+++ b/tutorials/multiBeam/system/blockMeshDict
@@ -35,7 +35,7 @@ vertices
 
 blocks
 (
-    hex (0 1 2 3 4 5 6 7) (300 50 30) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) (150 25 15) simpleGrading (1 1 1)
 );
 
 edges
diff --git a/tutorials/multiBeam/system/controlDict b/tutorials/multiBeam/system/controlDict
index 9657574..1cd9e90 100644
--- a/tutorials/multiBeam/system/controlDict
+++ b/tutorials/multiBeam/system/controlDict
@@ -7,11 +7,11 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      controlDict;
+    version         2;
+    format          ascii;
+    class           dictionary;
+    location        "system";
+    object          controlDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -23,9 +23,9 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.0014;
+endTime         0.0042;
 
-deltaT          1e-7;
+deltaT          1e-07;
 
 writeControl    adjustableRunTime;
 
@@ -48,6 +48,15 @@ runTimeModifiable yes;
 adjustTimeStep  yes;
 
 maxCo           0.5;
+
 maxDi           10;
+
 maxAlphaCo      1;
+
+functions
+{
+    #includeFunc    ExaCA
+}
+
+
 // ************************************************************************* //
diff --git a/tutorials/multiLayerPBF/Allclean b/tutorials/multiLayerPBF/Allclean
index 6b45f82..98bbd51 100755
--- a/tutorials/multiLayerPBF/Allclean
+++ b/tutorials/multiLayerPBF/Allclean
@@ -6,4 +6,9 @@ cd ${0%/*} || exit 1    # Run from this directory
 
 cleanCase
 
+rm -rf layer*
 #------------------------------------------------------------------------------
+rm -rf ExaCA/Output*
+rm -rf ExaCA/time-temperature.csv
+#------------------------------------------------------------------------------
+
diff --git a/tutorials/multiLayerPBF/Allrun b/tutorials/multiLayerPBF/Allrun
index 4fa5791..1503c12 100755
--- a/tutorials/multiLayerPBF/Allrun
+++ b/tutorials/multiLayerPBF/Allrun
@@ -1,28 +1,47 @@
 #!/bin/sh
-#------------------------------------------------------------------------------
-# =========                 |
-# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-#  \\    /   O peration     | Website:  https://openfoam.org
-#   \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
-#    \\/     M anipulation  |
-#------------------------------------------------------------------------------
-# Script
-#     Allrun for multiLayer
-#
-# Description
-#   Allrun script created for calling multiLayer utilities.
-#   User must specify path for multiLayer applications
-#
-#------------------------------------------------------------------------------
-cd ${0%/*} || exit 1
+cd ${0%/*} || exit 1    # Run from this directory
 
+# Source tutorial run functions
 . $WM_PROJECT_DIR/bin/tools/RunFunctions
+application=`getApplication`
+
+# Parse arguments
+withExaCA=false
+while [ "$#" -gt 0 ]; do
+  case "$1" in
+    -withExaCA)
+      withExaCA=true
+      ;;
+  esac
+  shift
+done
 
+if [ "$withExaCA" = true ]; then
+    foamDictionary -entry functions -set "{#includeFunc ExaCA}" system/controlDict
+else
+    foamDictionary -entry functions -set "{}" system/controlDict
+fi
+
+#------------------------------------------------------------------------------
+# AdditiveFOAM
 runApplication blockMesh
 
-path="$WM_PROJECT_INST_DIR/openfoam-am-dev/applications/utilities/multiLayer"
+# Change to AdditiveFOAM install path
+path="$WM_PROJECT_INST_DIR/AdditiveFOAM/applications/utilities/multiLayer"
 
-$path/runLayers -nLayers 5 -nCellsPerLayer 2 -layerThickness 20e-6
+$path/runLayers -nLayers 4 -nCellsPerLayer 4 -layerThickness 40e-6
 $path/reconstructLayers
 
 #------------------------------------------------------------------------------
+# ExaCA (Version 2.0.1)
+if [ "$withExaCA" = true ]; then
+    for d in layer*; do
+        echo "x,y,z,tm,ts,cr" > "$d/ExaCA/time-temperature.csv"
+        cat "$d/ExaCA/data_"* >> "$d/ExaCA/time-temperature.csv"
+        rm -rf "$d/ExaCA/data_*"
+    done
+    
+    NPROCS="$(foamDictionary -expand -entry numberOfSubdomains -value system/decomposeParDict)"
+    mpirun -np $NPROCS ~/install/exaca/bin/ExaCA ExaCA/input.json
+fi
+#------------------------------------------------------------------------------
diff --git a/tutorials/multiLayerPBF/ExaCA/input.json b/tutorials/multiLayerPBF/ExaCA/input.json
new file mode 100644
index 0000000..01de868
--- /dev/null
+++ b/tutorials/multiLayerPBF/ExaCA/input.json
@@ -0,0 +1,44 @@
+{
+   "SimulationType": "FromFile",
+   "MaterialFileName": "Inconel625.json",
+   "GrainOrientationFile": "GrainOrientationVectors.csv",
+   "RandomSeed": 0,
+   "Domain": {
+       "CellSize": 2.5,
+       "TimeStep": 0.125,
+       "NumberOfLayers": 4,
+       "LayerOffset": 16
+   },
+   "Nucleation": {
+      "Density": 100,
+      "MeanUndercooling": 5,
+      "StDev": 0.5
+   },
+   "TemperatureData": {
+       "LayerwiseTempRead": true,
+       "TemperatureFiles": [
+        "./layer0/ExaCA/time-temperature.csv",
+        "./layer1/ExaCA/time-temperature.csv",
+        "./layer2/ExaCA/time-temperature.csv",
+        "./layer3/ExaCA/time-temperature.csv"]
+   },
+   "Substrate": {
+      "MeanBaseplateGrainSize": 8.3,
+      "MeanPowderGrainSize": 2.5
+   },
+   "Printing": {
+      "PathToOutput": "./ExaCA/",
+      "OutputFile": "Output",
+      "PrintBinary": false,
+      "PrintExaConstitSize": 0,
+      "Interlayer": {
+          "Fields": [
+            "GrainID",
+            "LayerID",
+            "GrainMisorientation",
+            "UndercoolingCurrent",
+            "SolidificationEventCounter"],
+          "Layers": [0, 1]
+      }
+   }
+}
diff --git a/tutorials/multiLayerPBF/constant/heatSourceDict b/tutorials/multiLayerPBF/constant/heatSourceDict
index 0122d92..0c24aa0 100644
--- a/tutorials/multiLayerPBF/constant/heatSourceDict
+++ b/tutorials/multiLayerPBF/constant/heatSourceDict
@@ -17,23 +17,32 @@ FoamFile
 
 sources (beam);
 
+/*---------------------------------------------------------------------------*\
+Citation:
+    J. Coleman, G.L. Knapp, B. Stump, M. Rolchigo, K. Kincaid, A. Plotkowski,
+    A dynamic volumetric heat source model for laser additive manufacturing,
+    Additive Manufacturing (2024), https://doi.org/10.1016/j.addma.2024.104531.
+\*---------------------------------------------------------------------------*/
 beam
 {
-    pathName            scanPath;    
-
-    absorptionModel     constant;
-    
-    constantCoeffs
+    pathName            scanPath;
+ 
+    absorptionModel     Kelly;    
+    KellyCoeffs
     {
-        eta             0.3;
-    }
-    
-    heatSourceModel     superGaussian;
-    
-    superGaussianCoeffs
+        geometry        cone;
+        eta0            0.28;
+        etaMin          0.35;
+    }   
+
+    heatSourceModel     modifiedSuperGaussian;
+    modifiedSuperGaussianCoeffs
     {
-        k               2.0;
-        dimensions      (85.0e-6 85.0e-6 30e-6);
+        dimensions      (40.0e-6 40.0e-6 30e-6);
+        m               2.72;
+        k               7.95;
+        transient       true;
+        isoValue        1620;
         nPoints         (10 10 10);
     }
 }
diff --git a/tutorials/multiLayerPBF/constant/scanPath b/tutorials/multiLayerPBF/constant/scanPath
index b89b0ce..5fe3a9f 100644
--- a/tutorials/multiLayerPBF/constant/scanPath
+++ b/tutorials/multiLayerPBF/constant/scanPath
@@ -1,3 +1,3 @@
 Mode    X       Y       Z   Power   Param
 1       0.000   0.000   0   0       0
-0       0.002   0.000   0   179.2   0.8
+0       0.002   0.000   0   195.0   0.8
diff --git a/tutorials/multiLayerPBF/constant/transportProperties b/tutorials/multiLayerPBF/constant/transportProperties
index b945315..4ed515a 100644
--- a/tutorials/multiLayerPBF/constant/transportProperties
+++ b/tutorials/multiLayerPBF/constant/transportProperties
@@ -29,8 +29,8 @@ liquid
 
 powder
 {
-    kappa   (-0.07707   0.00075   0.0);
-    Cp      (747.568    0.0       0.0);
+    kappa   (0.21767  0.00102   0.0);
+    Cp      (317.53   0.0       0.0);
 }
 
 //- fluid flow properties
diff --git a/tutorials/multiLayerPBF/system/ExaCA b/tutorials/multiLayerPBF/system/ExaCA
new file mode 100644
index 0000000..fb9b4c4
--- /dev/null
+++ b/tutorials/multiLayerPBF/system/ExaCA
@@ -0,0 +1,19 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  10
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+
+ExaCA
+{
+    type ExaCA;
+    libs ("libExaCAFunctionObject.so");
+    
+    box         (0 -0.0001 -0.0002) (0.002 0.0001 0);
+    dx          2.5e-6;
+    isoValue    1620;
+}
+
+// *********************************************************************** //
diff --git a/tutorials/multiLayerPBF/system/blockMeshDict b/tutorials/multiLayerPBF/system/blockMeshDict
index c96b440..cd3a132 100644
--- a/tutorials/multiLayerPBF/system/blockMeshDict
+++ b/tutorials/multiLayerPBF/system/blockMeshDict
@@ -35,7 +35,7 @@ vertices
 
 blocks
 (
-    hex (0 1 2 3 4 5 6 7) (300 50 30) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) (150 25 30) simpleGrading (1 1 1)
 );
 
 edges
diff --git a/tutorials/multiLayerPBF/system/controlDict b/tutorials/multiLayerPBF/system/controlDict
index 6db0b85..7a3443c 100644
--- a/tutorials/multiLayerPBF/system/controlDict
+++ b/tutorials/multiLayerPBF/system/controlDict
@@ -7,11 +7,11 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      controlDict;
+    version         2;
+    format          ascii;
+    class           dictionary;
+    location        "system";
+    object          controlDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -25,7 +25,7 @@ stopAt          endTime;
 
 endTime         0.004;
 
-deltaT          1e-7;
+deltaT          1e-07;
 
 writeControl    adjustableRunTime;
 
@@ -48,6 +48,15 @@ runTimeModifiable yes;
 adjustTimeStep  yes;
 
 maxCo           0.5;
+
 maxDi           10;
+
 maxAlphaCo      1;
+
+functions
+{
+    #includeFunc    ExaCA
+}
+
+
 // ************************************************************************* //
diff --git a/tutorials/multiLayerPBF/system/decomposeParDict b/tutorials/multiLayerPBF/system/decomposeParDict
index fe538dd..08d37ba 100644
--- a/tutorials/multiLayerPBF/system/decomposeParDict
+++ b/tutorials/multiLayerPBF/system/decomposeParDict
@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains  4;
+numberOfSubdomains  6;
 
 method              scotch;
 
diff --git a/tutorials/multiLayerPBF/system/setFieldsDict b/tutorials/multiLayerPBF/system/setFieldsDict
index 7cfaa06..601f94b 100644
--- a/tutorials/multiLayerPBF/system/setFieldsDict
+++ b/tutorials/multiLayerPBF/system/setFieldsDict
@@ -7,11 +7,11 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      setFieldsDict;
+    version         2;
+    format          ascii;
+    class           dictionary;
+    location        "system";
+    object          setFieldsDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -19,7 +19,7 @@ regions
 (
     boxToCell
     {
-        box (-1 -1 -2e-5) (1 1 1);
+        box (-1 -1 -4e-5) (1 1 1);
 
         fieldValues
         (
diff --git a/tutorials/marangoni/0/T b/tutorials/verificationAndValidation/marangoni/0/T
similarity index 100%
rename from tutorials/marangoni/0/T
rename to tutorials/verificationAndValidation/marangoni/0/T
diff --git a/tutorials/marangoni/0/U b/tutorials/verificationAndValidation/marangoni/0/U
similarity index 100%
rename from tutorials/marangoni/0/U
rename to tutorials/verificationAndValidation/marangoni/0/U
diff --git a/tutorials/marangoni/0/alpha.solid b/tutorials/verificationAndValidation/marangoni/0/alpha.solid
similarity index 100%
rename from tutorials/marangoni/0/alpha.solid
rename to tutorials/verificationAndValidation/marangoni/0/alpha.solid
diff --git a/tutorials/marangoni/0/p_rgh b/tutorials/verificationAndValidation/marangoni/0/p_rgh
similarity index 100%
rename from tutorials/marangoni/0/p_rgh
rename to tutorials/verificationAndValidation/marangoni/0/p_rgh
diff --git a/tutorials/marangoni/Allclean b/tutorials/verificationAndValidation/marangoni/Allclean
similarity index 100%
rename from tutorials/marangoni/Allclean
rename to tutorials/verificationAndValidation/marangoni/Allclean
diff --git a/tutorials/marangoni/Allrun b/tutorials/verificationAndValidation/marangoni/Allrun
similarity index 100%
rename from tutorials/marangoni/Allrun
rename to tutorials/verificationAndValidation/marangoni/Allrun
diff --git a/tutorials/marangoni/constant/g b/tutorials/verificationAndValidation/marangoni/constant/g
similarity index 100%
rename from tutorials/marangoni/constant/g
rename to tutorials/verificationAndValidation/marangoni/constant/g
diff --git a/tutorials/marangoni/constant/heatSourceDict b/tutorials/verificationAndValidation/marangoni/constant/heatSourceDict
similarity index 100%
rename from tutorials/marangoni/constant/heatSourceDict
rename to tutorials/verificationAndValidation/marangoni/constant/heatSourceDict
diff --git a/tutorials/marangoni/constant/scanPath b/tutorials/verificationAndValidation/marangoni/constant/scanPath
similarity index 100%
rename from tutorials/marangoni/constant/scanPath
rename to tutorials/verificationAndValidation/marangoni/constant/scanPath
diff --git a/tutorials/marangoni/constant/thermoPath b/tutorials/verificationAndValidation/marangoni/constant/thermoPath
similarity index 100%
rename from tutorials/marangoni/constant/thermoPath
rename to tutorials/verificationAndValidation/marangoni/constant/thermoPath
diff --git a/tutorials/marangoni/constant/transportProperties b/tutorials/verificationAndValidation/marangoni/constant/transportProperties
similarity index 100%
rename from tutorials/marangoni/constant/transportProperties
rename to tutorials/verificationAndValidation/marangoni/constant/transportProperties
diff --git a/tutorials/marangoni/system/blockMeshDict b/tutorials/verificationAndValidation/marangoni/system/blockMeshDict
similarity index 100%
rename from tutorials/marangoni/system/blockMeshDict
rename to tutorials/verificationAndValidation/marangoni/system/blockMeshDict
diff --git a/tutorials/marangoni/system/controlDict b/tutorials/verificationAndValidation/marangoni/system/controlDict
similarity index 100%
rename from tutorials/marangoni/system/controlDict
rename to tutorials/verificationAndValidation/marangoni/system/controlDict
diff --git a/tutorials/marangoni/system/decomposeParDict b/tutorials/verificationAndValidation/marangoni/system/decomposeParDict
similarity index 100%
rename from tutorials/marangoni/system/decomposeParDict
rename to tutorials/verificationAndValidation/marangoni/system/decomposeParDict
diff --git a/tutorials/marangoni/system/fvSchemes b/tutorials/verificationAndValidation/marangoni/system/fvSchemes
similarity index 100%
rename from tutorials/marangoni/system/fvSchemes
rename to tutorials/verificationAndValidation/marangoni/system/fvSchemes
diff --git a/tutorials/marangoni/system/fvSolution b/tutorials/verificationAndValidation/marangoni/system/fvSolution
similarity index 100%
rename from tutorials/marangoni/system/fvSolution
rename to tutorials/verificationAndValidation/marangoni/system/fvSolution