diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py
index 0f1c2873..a83cfb96 100755
--- a/dafoam/pyDAFoam.py
+++ b/dafoam/pyDAFoam.py
@@ -508,6 +508,17 @@ def __init__(self):
         ## constrainHbyA back to the primal and adjoint solvers.
         self.useConstrainHbyA = False
 
+        ## parameters for regression models
+        self.regressionModel = {
+            "active": False,
+            "modelType": "neuralNetwork",
+            "inputNames": ["None"],
+            "outputName": "None",
+            "hiddenLayerNeurons": [0],
+            "outputShift": 0.0,
+            "outputScale": 1.0,
+        }
+
         # *********************************************************************************************
         # ************************************ Advance Options ****************************************
         # *********************************************************************************************
@@ -2264,6 +2275,44 @@ def calcTotalDerivsACT(self, objFuncName, designVarName, designVarType, dFScalin
             else:
                 self.adjTotalDeriv[objFuncName][designVarName][i] = totalDerivSeq[i]
 
+    def calcTotalDerivsRegPar(self, objFuncName, designVarName, accumulateTotal=False):
+
+        xDV = self.DVGeo.getValues()
+        nDVs = len(xDV[designVarName])
+
+        nParameters = self.solver.getNRegressionParameters()
+        if nDVs != nParameters:
+            raise Error("number of parameters not valid!")
+
+        # We assume  dFdRegPar is always zero
+        # call the total deriv
+
+        xvArray = self.vec2Array(self.xvVec)
+        wArray = self.vec2Array(self.wVec)
+        seedArray = self.vec2Array(self.adjVectors[objFuncName])
+        parameters = xDV[designVarName].copy(order="C")
+        productArray = np.zeros(nDVs)
+
+        # calculate dRdFieldT*Psi and save it to totalDeriv
+        self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, productArray)
+        # all reduce because parameters is a global DV
+        productArray = self.comm.allreduce(productArray, op=MPI.SUM)
+
+        # totalDeriv = dFdRegPar - dRdRegParT*psi
+        productArray *= -1.0
+
+        # assign the total derivative to self.adjTotalDeriv
+        if self.adjTotalDeriv[objFuncName][designVarName] is None:
+            self.adjTotalDeriv[objFuncName][designVarName] = np.zeros(nDVs, self.dtype)
+
+        # NOTE: productArray is already in Seq
+
+        for i in range(nDVs):
+            if accumulateTotal is True:
+                self.adjTotalDeriv[objFuncName][designVarName][i] += productArray[i]
+            else:
+                self.adjTotalDeriv[objFuncName][designVarName][i] = productArray[i]
+
     def solveAdjointUnsteady(self):
         """
         Run adjoint solver to compute the total derivs for unsteady solvers
@@ -2274,6 +2323,8 @@ def solveAdjointUnsteady(self):
 
         wVec: vector that contains all the state variables
 
+        designVariable: a dictionary that contain the design variable values
+
         Output:
         -------
         self.adjTotalDeriv: the dict contains all the total derivative vectors
@@ -2462,6 +2513,8 @@ def solveAdjointUnsteady(self):
                         elif designVarDict[designVarName]["designVarType"] == "Field":
                             fieldType = designVarDict[designVarName]["fieldType"]
                             self.calcTotalDerivsField(objFuncName, designVarName, fieldType, dFScaling, True)
+                        elif designVarDict[designVarName]["designVarType"] == "RegPar":
+                            self.calcTotalDerivsRegPar(objFuncName, designVarName, True)
                         else:
                             raise Error("designVarType not valid!")
 
@@ -2751,6 +2804,15 @@ def solveAdjoint(self):
                             self.calcTotalDerivsField(objFuncName, designVarName, fieldType)
                 else:
                     raise Error("For Field design variable type, we only support useAD->mode=reverse")
+            ################### RegPar: regression model's parameters ###################
+            elif designVarDict[designVarName]["designVarType"] == "RegPar":
+                if self.getOption("useAD")["mode"] == "reverse":
+                    # loop over all objectives
+                    for objFuncName in objFuncDict:
+                        if objFuncName in self.objFuncNames4Adj:
+                            self.calcTotalDerivsRegPar(objFuncName, designVarName)
+                else:
+                    raise Error("For Field design variable type, we only support useAD->mode=reverse")
             else:
                 raise Error("designVarType %s not supported!" % designVarDict[designVarName]["designVarType"])
 
@@ -3151,7 +3213,7 @@ def _initSolver(self):
         adjMode = self.getOption("unsteadyAdjoint")["mode"]
         if adjMode == "hybrid":
             self.initTimeInstanceMats()
-        
+
         Info("Init solver done! ElapsedClockTime %f s" % self.solver.getElapsedClockTime())
         Info("Init solver done!. ElapsedCpuTime %f s" % self.solver.getElapsedCpuTime())
 
@@ -4002,6 +4064,23 @@ def _initOption(self, name, value):
                 "Expected data type is %-47s \n "
                 "Received data type is %-47s" % (name, self.defaultOptions[name][0], type(value))
             )
+ 
+    def setRegressionParameter(self, idx, val):
+        """
+        Update the regression parameters
+
+        Parameters
+        ----------
+        idx: int
+            the index of the parameter
+
+        val, float
+            the parameter value to set
+        """
+
+        self.solver.setRegressionParameter(idx, val)
+        if self.getOption("useAD")["mode"] in ["forward", "reverse"]:
+            self.solverAD.setRegressionParameter(idx, val)
 
     def setFieldValue4GlobalCellI(self, fieldName, val, globalCellI, compI=0):
         """
diff --git a/src/adjoint/DAModel/DAModel.C b/src/adjoint/DAModel/DAModel.C
index 4fe38cf1..91c93cd9 100755
--- a/src/adjoint/DAModel/DAModel.C
+++ b/src/adjoint/DAModel/DAModel.C
@@ -292,6 +292,28 @@ void DAModel::getTurbProdTerm(scalarList& prodTerm) const
 
 }
 
+void DAModel::getTurbProdOverDestruct(scalarList& PoD) const
+{
+    /*
+    Description: 
+        Return the value of the production/destruction term from the turbulence model 
+    */
+
+#ifndef SolidDASolver
+    if (hasTurbulenceModel_)
+    {
+        DATurbulenceModel& daTurb = const_cast<DATurbulenceModel&>(
+            mesh_.thisDb().lookupObject<DATurbulenceModel>("DATurbulenceModel"));
+        daTurb.getTurbProdOverDestruct(PoD);
+    }
+
+    if (hasRadiationModel_)
+    {
+    }
+#endif
+
+}
+
 #ifdef CompressibleFlow
 const fluidThermo& DAModel::getThermo() const
 {
diff --git a/src/adjoint/DAModel/DAModel.H b/src/adjoint/DAModel/DAModel.H
index 5c074378..ba2b610e 100755
--- a/src/adjoint/DAModel/DAModel.H
+++ b/src/adjoint/DAModel/DAModel.H
@@ -105,9 +105,12 @@ public:
     /// update intermediate variables that are dependent on the model states
     void updateIntermediateVariables();
 
-    /// return the value of the production term from the turbulence model 
+    /// return the value of the production term from the turbulence model
     void getTurbProdTerm(scalarList& prodTerm) const;
 
+    /// return the value of the destruction term from the turbulence model
+    void getTurbProdOverDestruct(scalarList& PoD) const;
+
 #ifndef SolidDASolver
     /// get a reference to DATurbulenceModel
     const DATurbulenceModel& getDATurbulenceModel() const;
@@ -117,7 +120,6 @@ public:
     /// return thermo_ from DATurbulenceModel
     const fluidThermo& getThermo() const;
 #endif
-
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C
index 977ba162..16de8cec 100755
--- a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C
+++ b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C
@@ -518,6 +518,27 @@ void DASpalartAllmaras::getFvMatrixFields(
     upper = nuTildaEqn.upper();
     lower = nuTildaEqn.lower();
 }
+
+void DASpalartAllmaras::getTurbProdOverDestruct(scalarList& PoD) const
+{
+    /*
+    Description:
+        Return the value of the production over destruction term from the turbulence model 
+    */
+
+    const volScalarField chi(this->chi());
+    const volScalarField fv1(this->fv1(chi));
+
+    const volScalarField Stilda(this->Stilda(chi, fv1));
+
+    volScalarField P = Cb1_ * phase_ * rho_ * Stilda * nuTilda_;
+    volScalarField D = Cw1_ * phase_ * rho_ * fw(Stilda) * sqr(nuTilda_ / y_);
+
+    forAll(P, cellI)
+    {
+        PoD[cellI] = P[cellI] / D[cellI];
+    }
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H
index befc48f1..1f1b485a 100755
--- a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H
+++ b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H
@@ -146,6 +146,9 @@ public:
         scalarField& diag,
         scalarField& upper,
         scalarField& lower);
+    
+    /// return the value of the destruction term from the turbulence model
+    virtual void getTurbProdOverDestruct(scalarList& PoD) const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C
index 97a6edc9..8382c61e 100755
--- a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C
+++ b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C
@@ -737,13 +737,36 @@ void DASpalartAllmarasFv3::getFvMatrixFields(
             - Cb2_ / sigmaNut_ * phase_ * rho_ * magSqr(fvc::grad(nuTilda_))
         == Cb1_ * phase_ * rho_ * Stilda * nuTilda_ * betaFI_
             - fvm::Sp(Cw1_ * phase_ * rho_ * fw(Stilda) * nuTilda_ / sqr(y_), nuTilda_));
-    
+
     nuTildaEqn.relax();
 
     diag = nuTildaEqn.D();
     upper = nuTildaEqn.upper();
     lower = nuTildaEqn.lower();
 }
+
+void DASpalartAllmarasFv3::getTurbProdOverDestruct(scalarList& PoD) const
+{
+    /*
+    Description:
+        Return the value of the production over destruction term from the turbulence model 
+    */
+
+    const volScalarField chi(this->chi());
+    const volScalarField fv1(this->fv1(chi));
+
+    const volScalarField Stilda(
+        this->fv3(chi, fv1) * ::sqrt(2.0) * mag(skew(fvc::grad(U_)))
+        + this->fv2(chi, fv1) * nuTilda_ / sqr(kappa_ * y_));
+
+    volScalarField P = Cb1_ * phase_ * rho_ * Stilda * nuTilda_;
+    volScalarField D = Cw1_ * phase_ * rho_ * fw(Stilda) * sqr(nuTilda_ / y_);
+
+    forAll(P, cellI)
+    {
+        PoD[cellI] = P[cellI] / D[cellI];
+    }
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H
index 91d2fef4..09702d98 100755
--- a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H
+++ b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H
@@ -168,6 +168,9 @@ public:
         scalarField& diag,
         scalarField& upper,
         scalarField& lower);
+
+    /// return the value of the destruction term from the turbulence model
+    virtual void getTurbProdOverDestruct(scalarList& PoD) const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C b/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C
index 2ed04794..71f81b48 100755
--- a/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C
+++ b/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C
@@ -487,7 +487,19 @@ void DATurbulenceModel::getTurbProdTerm(scalarList& prodTerm) const
         Return the value of the production term from the turbulence model 
     */
 
-    FatalErrorIn("DATurbulenceModel::getSAProdTerm")
+    FatalErrorIn("DATurbulenceModel::getTurbProdTerm")
+        << "Child class not implemented!"
+        << abort(FatalError);
+}
+
+void DATurbulenceModel::getTurbProdOverDestruct(scalarList& PoD) const
+{
+    /*
+    Description:
+        Return the value of the production over destruction term from the turbulence model 
+    */
+
+    FatalErrorIn("DATurbulenceModel::getTurbProdOverDestruct")
         << "Child class not implemented!"
         << abort(FatalError);
 }
diff --git a/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H b/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H
index dc936934..e08bb15c 100755
--- a/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H
+++ b/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H
@@ -211,9 +211,12 @@ public:
     /// solve the residual equations and update the state
     virtual void correct(label printToScreen) = 0;
 
-    /// return the value of the production term from the Spalart Allmaras model
+    /// return the value of the production term from the turbulence model
     virtual void getTurbProdTerm(scalarList& prodTerm) const;
 
+    /// return the value of the destruction term from the turbulence model
+    virtual void getTurbProdOverDestruct(scalarList& PoD) const;
+
     /// dev terms
     tmp<volSymmTensorField> devRhoReff() const;
 
diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C b/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C
index 0ec19b99..3c9b3718 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C
@@ -51,6 +51,8 @@ DAObjFuncVariableVolSum::DAObjFuncVariableVolSum(
 
     objFuncDict_.readEntry<label>("divByTotalVol", divByTotalVol_);
 
+    multiplyVol_ = objFuncDict_.lookupOrDefault<label>("multiplyVol", 1);
+
     if (daIndex.adjStateNames.found(varName_))
     {
         objFuncConInfo_ = {{varName_}};
@@ -114,7 +116,11 @@ void DAObjFuncVariableVolSum::calcObjFunc(
         forAll(objFuncCellSources, idxI)
         {
             const label& cellI = objFuncCellSources[idxI];
-            scalar volume = mesh_.V()[cellI];
+            scalar volume = 1.0;
+            if (multiplyVol_)
+            {
+                volume = mesh_.V()[cellI];
+            }
             if (isSquare_)
             {
                 objFuncCellValues[idxI] = scale_ * volume * var[cellI] * var[cellI];
@@ -133,7 +139,11 @@ void DAObjFuncVariableVolSum::calcObjFunc(
         forAll(objFuncCellSources, idxI)
         {
             const label& cellI = objFuncCellSources[idxI];
-            scalar volume = mesh_.V()[cellI];
+            scalar volume = 1.0;
+            if (multiplyVol_)
+            {
+                volume = mesh_.V()[cellI];
+            }
             if (isSquare_)
             {
                 objFuncCellValues[idxI] = scale_ * volume * var[cellI][component_] * var[cellI][component_];
diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.H b/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.H
index a6e07456..e4a3d25f 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.H
+++ b/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.H
@@ -44,6 +44,9 @@ protected:
     /// whether to normalize the volume sum with the total volume
     label divByTotalVol_;
 
+    /// whether to multiply the variable by the volume
+    label multiplyVol_;
+
 public:
     TypeName("variableVolSum");
     // Constructors
diff --git a/src/adjoint/DARegression/DARegression.C b/src/adjoint/DARegression/DARegression.C
new file mode 100755
index 00000000..5c9c1308
--- /dev/null
+++ b/src/adjoint/DARegression/DARegression.C
@@ -0,0 +1,253 @@
+/*---------------------------------------------------------------------------*\
+
+    DAFoam  : Discrete Adjoint with OpenFOAM
+    Version : v3
+
+\*---------------------------------------------------------------------------*/
+
+#include "DARegression.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+DARegression::DARegression(
+    const fvMesh& mesh,
+    const DAOption& daOption,
+    const DAModel& daModel)
+    : mesh_(mesh),
+      daOption_(daOption),
+      daModel_(daModel)
+{
+    dictionary regSubDict = daOption.getAllOptions().subDict("regressionModel");
+    regSubDict.readEntry<word>("modelType", modelType_);
+    regSubDict.readEntry<wordList>("inputNames", inputNames_);
+    regSubDict.readEntry<word>("outputName", outputName_);
+
+    regSubDict.readEntry<scalar>("outputShift", outputShift_);
+
+    regSubDict.readEntry<scalar>("outputScale", outputScale_);
+
+    active_ = regSubDict.getLabel("active");
+
+    if (modelType_ == "neuralNetwork")
+    {
+        regSubDict.readEntry<labelList>("hiddenLayerNeurons", hiddenLayerNeurons_);
+    }
+
+    // initialize parameters and give it large values
+    if (active_)
+    {
+        label nParameters = this->nParameters();
+        parameters_.setSize(nParameters);
+        forAll(parameters_, idxI)
+        {
+            parameters_[idxI] = 1e16;
+        }
+    }
+}
+
+void DARegression::compute()
+{
+    /*
+    Description:
+        Calculate the prescribed output field based on the prescribed input fields using a regression model.
+        We support only neural network at this moment
+
+    Input:
+        parameters_:
+            the parameters for the regression. For the neural network model, these parameters
+            are the weights and biases. We order them in an 1D array, starting from the first input's weight and biases.
+
+        daOption Dict: regressionModel
+            inputNames: a list of volScalarFields prescribed by inputFields
+            hiddenLayerNeurons: number of neurons for each hidden layer. 
+            example: {5, 3, 4} means three hidden layers with 5, 3, and 4 neurons.
+
+    Output:
+        a volScalarField prescribed by outputName
+    */
+
+    if (!active_)
+    {
+        return;
+    }
+
+    volScalarField& outputField = const_cast<volScalarField&>(mesh_.thisDb().lookupObject<volScalarField>(outputName_));
+
+    if (modelType_ == "neuralNetwork")
+    {
+        label nHiddenLayers = hiddenLayerNeurons_.size();
+
+        List<List<scalar>> inputFields;
+        inputFields.setSize(inputNames_.size());
+
+        forAll(inputNames_, idxI)
+        {
+            inputFields[idxI].setSize(mesh_.nCells());
+            word inputName = inputNames_[idxI];
+            if (inputName == "SoQ")
+            {
+                // Shear / vorticity
+                const volVectorField& U = mesh_.thisDb().lookupObject<volVectorField>("U");
+                const tmp<volTensorField> tgradU(fvc::grad(U));
+                const volTensorField& gradU = tgradU();
+                volScalarField magOmega = mag(skew(gradU));
+                volScalarField magS = mag(symm(gradU));
+                forAll(inputFields[idxI], cellI)
+                {
+                    inputFields[idxI][cellI] = magS[cellI] / magOmega[cellI];
+                }
+            }
+            else if (inputName == "PoD")
+            {
+                // production / destruction
+                daModel_.getTurbProdOverDestruct(inputFields[idxI]);
+            }
+            else if (inputName == "chiSA")
+            {
+#ifndef SolidDASolver
+                // the chi() function from SA
+                const volScalarField& nuTilda = mesh_.thisDb().lookupObject<volScalarField>("nuTilda");
+                volScalarField nu = daModel_.getDATurbulenceModel().nu();
+                forAll(inputFields[idxI], cellI)
+                {
+                    inputFields[idxI][cellI] = nuTilda[cellI] / nu[cellI];
+                }
+#endif
+            }
+            else if (inputName == "pGradStream")
+            {
+                // pressure gradient along stream
+                const volScalarField& p = mesh_.thisDb().lookupObject<volScalarField>("p");
+                const volVectorField& U = mesh_.thisDb().lookupObject<volVectorField>("U");
+                volVectorField pGrad("gradP", fvc::grad(p));
+                volScalarField pG_denominator(mag(U) * mag(pGrad) + mag(U & pGrad));
+                volScalarField pGradAlongStream = (U & pGrad) / Foam::max(pG_denominator, dimensionedScalar("minpG", dimensionSet(0, 2, -3, 0, 0, 0, 0), SMALL));
+                forAll(inputFields[idxI], cellI)
+                {
+                    inputFields[idxI][cellI] = pGradAlongStream[cellI];
+                }
+            }
+            else
+            {
+                FatalErrorIn("") << "inputName: " << inputName << " not supported. Options are: SoQ, PoD, chiSA, pGradStream" << abort(FatalError);
+            }
+        }
+
+        List<List<scalar>> layerVals;
+        layerVals.setSize(nHiddenLayers);
+        for (label layerI = 0; layerI < nHiddenLayers; layerI++)
+        {
+            label nNeurons = hiddenLayerNeurons_[layerI];
+            layerVals[layerI].setSize(nNeurons);
+        }
+
+        forAll(mesh_.cells(), cellI)
+        {
+            label counterI = 0;
+
+            for (label layerI = 0; layerI < nHiddenLayers; layerI++)
+            {
+                label nNeurons = hiddenLayerNeurons_[layerI];
+                forAll(layerVals[layerI], neuronI)
+                {
+                    layerVals[layerI][neuronI] = 0.0;
+                }
+                for (label neuronI = 0; neuronI < nNeurons; neuronI++)
+                {
+                    if (layerI == 0)
+                    {
+                        forAll(inputNames_, neuronJ)
+                        {
+                            // weighted sum
+                            layerVals[layerI][neuronI] += inputFields[neuronJ][cellI] * parameters_[counterI];
+                            counterI++;
+                        }
+                    }
+                    else
+                    {
+                        forAll(layerVals[layerI - 1], neuronJ)
+                        {
+                            // weighted sum
+                            layerVals[layerI][neuronI] += layerVals[layerI - 1][neuronJ] * parameters_[counterI];
+                            counterI++;
+                        }
+                    }
+                    // bias
+                    layerVals[layerI][neuronI] += parameters_[counterI];
+                    counterI++;
+                    // activation function
+                    layerVals[layerI][neuronI] = 1 / (1 + exp(-layerVals[layerI][neuronI]));
+                }
+            }
+            // final output layer, we have only one output
+            scalar outputVal = 0.0;
+            forAll(layerVals[nHiddenLayers - 1], neuronJ)
+            {
+                // weighted sum
+                outputVal += layerVals[nHiddenLayers - 1][neuronJ] * parameters_[counterI];
+                counterI++;
+            }
+            // bias
+            outputVal += parameters_[counterI];
+
+            outputField[cellI] = outputScale_ * (outputVal + outputShift_);
+        }
+
+        outputField.correctBoundaryConditions();
+    }
+    else
+    {
+        FatalErrorIn("") << "modelType_: " << modelType_ << " not supported. Options are: neuralNetwork" << abort(FatalError);
+    }
+}
+
+label DARegression::nParameters()
+{
+    /*
+    Description:
+        get the number of parameters
+    */
+
+    if (!active_)
+    {
+        FatalErrorIn("") << "nParameters() is called but the regression model is not active!" << abort(FatalError);
+    }
+
+    if (modelType_ == "neuralNetwork")
+    {
+        label nHiddenLayers = hiddenLayerNeurons_.size();
+        label nInputs = inputNames_.size();
+
+        // add weights
+        // input
+        label nParameters = nInputs * hiddenLayerNeurons_[0];
+        // hidden layers
+        for (label layerI = 1; layerI < hiddenLayerNeurons_.size(); layerI++)
+        {
+            nParameters += hiddenLayerNeurons_[layerI] * hiddenLayerNeurons_[layerI - 1];
+        }
+        // output
+        nParameters += hiddenLayerNeurons_[nHiddenLayers - 1] * 1;
+
+        // add biases
+        // add hidden layers
+        for (label layerI = 0; layerI < hiddenLayerNeurons_.size(); layerI++)
+        {
+            nParameters += hiddenLayerNeurons_[layerI];
+        }
+        // add output layer
+        nParameters += 1;
+
+        return nParameters;
+    }
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/adjoint/DARegression/DARegression.H b/src/adjoint/DARegression/DARegression.H
new file mode 100755
index 00000000..9ff93e3a
--- /dev/null
+++ b/src/adjoint/DARegression/DARegression.H
@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+
+    DAFoam  : Discrete Adjoint with OpenFOAM
+    Version : v3
+
+    Description:
+        Regression model
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DARegression_H
+#define DARegression_H
+
+#include "fvOptions.H"
+#include "surfaceFields.H"
+#include "DAOption.H"
+#include "DAUtility.H"
+#include "DAModel.H"
+#include "globalIndex.H"
+#include "DAMacroFunctions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class DARegression Declaration
+\*---------------------------------------------------------------------------*/
+
+class DARegression
+{
+
+private:
+    /// Disallow default bitwise copy construct
+    DARegression(const DARegression&);
+
+    /// Disallow default bitwise assignment
+    void operator=(const DARegression&);
+
+protected:
+    /// Foam::fvMesh object
+    const fvMesh& mesh_;
+
+    /// Foam::DAOption object
+    const DAOption& daOption_;
+
+    /// DAModel object
+    const DAModel& daModel_;
+
+    /// the type of regression model
+    word modelType_;
+
+    /// a list of words for the inputs
+    wordList inputNames_;
+
+    /// a list of words for the outputs
+    word outputName_;
+
+    /// number of neurons hidden layers of the neural network
+    labelList hiddenLayerNeurons_;
+
+    /// we can shift the output. we always shift before scaling it.
+    scalar outputShift_;
+
+    /// we can scale the output. we always shift before scaling it.
+    scalar outputScale_;
+
+    /// the parameters for the regression model
+    scalarList parameters_;
+
+    /// whether the regression model is active
+    label active_;
+
+public:
+    /// Constructors
+    DARegression(
+        const fvMesh& mesh,
+        const DAOption& daOption,
+        const DAModel& daModel);
+
+    /// Destructor
+    virtual ~DARegression()
+    {
+    }
+
+    // Members
+
+    /// compute the output based on the latest parameters and inputs
+    void compute();
+
+    /// get the number of parameters for this regression model
+    label nParameters();
+
+    /// get a specific parameter value
+    scalar getParameter(label idxI)
+    {
+        return parameters_[idxI];
+    }
+
+    /// set a value to a parameter give an index and a value
+    void setParameter(label idxI, scalar val)
+    {
+        parameters_[idxI] = val;
+    }
+
+    /// check if the regression model is initialized properly
+    void validate()
+    {
+        if (!active_)
+        {
+            return;
+        }
+        else
+        {
+            scalar testVal = this->getParameter(0);
+            if (fabs(testVal) > 1e15)
+            {
+                FatalErrorIn("") << "regressionModel is active but the parameter values are not set!" << abort(FatalError);
+            }
+        }
+    }
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C
index 99331025..b1064886 100755
--- a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C
+++ b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C
@@ -175,6 +175,9 @@ label DAPimpleFoam::solvePrimal(
     scalar deltaT = runTime.deltaT().value();
     label nInstances = round(endTime / deltaT);
 
+    // check if the parameters are set in the Python layer
+    daRegressionPtr_->validate();
+
     // main loop
     for (label iter = 1; iter <= nInstances; iter++)
     {
@@ -208,6 +211,9 @@ label DAPimpleFoam::solvePrimal(
 #include "pEqnPimple.H"
             }
 
+            // update the output field value at each iteration, if the regression model is active
+            daRegressionPtr_->compute();
+
             laminarTransport.correct();
             daTurbulenceModelPtr_->correct(pimplePrintToScreen);
         }
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
index a1945f86..b681d2f1 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
@@ -146,6 +146,9 @@ label DARhoSimpleFoam::solvePrimal(
         }
     }
 
+    // check if the parameters are set in the Python layer
+    daRegressionPtr_->validate();
+
     primalMinRes_ = 1e10;
     label printInterval = daOptionPtr_->getOption<label>("printInterval");
     label printToScreen = 0;
@@ -167,6 +170,9 @@ label DARhoSimpleFoam::solvePrimal(
 #include "EEqnRhoSimple.H"
 #include "pEqnRhoSimple.H"
 
+        // update the output field value at each iteration, if the regression model is active
+        daRegressionPtr_->compute();
+
         daTurbulenceModelPtr_->correct(printToScreen);
 
         if (printToScreen)
diff --git a/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C b/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
index 5d52321b..ea3ca61b 100755
--- a/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
+++ b/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
@@ -102,6 +102,7 @@ void DASimpleFoam::initSolver()
         daFvSourcePtr_.reset(DAFvSource::New(
             fvSourceType, mesh, daOptionPtr_(), daModelPtr_(), daIndexPtr_()));
     }
+
 }
 
 label DASimpleFoam::solvePrimal(
@@ -156,6 +157,9 @@ label DASimpleFoam::solvePrimal(
         }
     }
 
+    // check if the parameters are set in the Python layer
+    daRegressionPtr_->validate();
+
     primalMinRes_ = 1e10;
     label printInterval = daOptionPtr_->getOption<label>("printInterval");
     label printToScreen = 0;
@@ -177,6 +181,9 @@ label DASimpleFoam::solvePrimal(
 #include "pEqnSimple.H"
         }
 
+        // update the output field value at each iteration, if the regression model is active
+        daRegressionPtr_->compute();
+
         laminarTransport.correct();
         daTurbulenceModelPtr_->correct(printToScreen);
 
@@ -1057,7 +1064,7 @@ void DASimpleFoam::invTranProdUEqn(
         We won't ADR this function, so we can treat most of the arguments as const
     */
 
-/*
+    /*
     const objectRegistry& db = meshPtr_->thisDb();
     const surfaceScalarField& phi = db.lookupObject<surfaceScalarField>("phi");
     volScalarField nuEff = daTurbulenceModelPtr_->nuEff();
@@ -1113,7 +1120,7 @@ void DASimpleFoam::invTranProdPEqn(
         We won't ADR this function, so we can treat most of the arguments as const
     */
 
-/*
+    /*
     const objectRegistry& db = meshPtr_->thisDb();
     const volVectorField& U = db.lookupObject<volVectorField>("U");
     const surfaceScalarField& phi = db.lookupObject<surfaceScalarField>("phi");
diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
index 57a9290c..5853a2cc 100644
--- a/src/adjoint/DASolver/DASolver.C
+++ b/src/adjoint/DASolver/DASolver.C
@@ -43,7 +43,8 @@ DASolver::DASolver(
       daFieldPtr_(nullptr),
       daCheckMeshPtr_(nullptr),
       daLinearEqnPtr_(nullptr),
-      daResidualPtr_(nullptr)
+      daResidualPtr_(nullptr),
+      daRegressionPtr_(nullptr)
 #ifdef CODI_AD_REVERSE
       ,
       globalADTape_(codi::RealReverse::getTape())
@@ -5866,6 +5867,125 @@ void DASolver::calcdRdThermalTPsiAD(
 #endif
 }
 
+void DASolver::calcdRdRegParTPsiAD(
+    const double* volCoords,
+    const double* states,
+    const double* parameters,
+    const double* seeds,
+    double* product)
+{
+#ifdef CODI_AD_REVERSE
+    /*
+    Description:
+        Compute the matrix-vector products [dR/dRegParameters]^T*Psi using reverse-mode AD
+    
+    Input:
+
+        xvVec: the volume mesh coordinate vector
+
+        wVec: the state variable vector
+
+        parameters: the parameters for the regression model
+
+        psi: the vector to multiply dRdXv
+    
+    Output:
+        prodVec: the matrix-vector products [dR/dRegParameters]^T*Psi
+    */
+
+    Info << "Calculating [dRdRegPar]^T * Psi using reverse-mode AD" << endl;
+
+    scalar* volCoordsArray = new scalar[daIndexPtr_->nLocalXv];
+    for (label i = 0; i < daIndexPtr_->nLocalXv; i++)
+    {
+        volCoordsArray[i] = volCoords[i];
+    }
+
+    scalar* statesArray = new scalar[daIndexPtr_->nLocalAdjointStates];
+    for (label i = 0; i < daIndexPtr_->nLocalAdjointStates; i++)
+    {
+        statesArray[i] = states[i];
+    }
+
+    label nParameters = this->getNRegressionParameters();
+    scalar* parametersArray = new scalar[nParameters];
+    for (label i = 0; i < nParameters; i++)
+    {
+        parametersArray[i] = parameters[i];
+    }
+
+    this->updateOFMesh(volCoordsArray);
+    this->updateOFField(statesArray);
+
+    // update the OpenFOAM variables and reset their seeds (gradient part) to zeros
+    this->resetOFSeeds();
+    // reset the AD tape
+    this->globalADTape_.reset();
+    // start recording
+    this->globalADTape_.setActive();
+
+    // register inputs
+    for (label i = 0; i < nParameters; i++)
+    {
+        this->globalADTape_.registerInput(parametersArray[i]);
+    }
+    for (label i = 0; i < nParameters; i++)
+    {
+        this->setRegressionParameter(i, parametersArray[i]);
+    }
+
+    // calculate outputs
+    this->regressionModelCompute();
+
+    // compute residuals
+    this->updateStateBoundaryConditions();
+    this->calcResiduals();
+
+    // register outputs
+    this->registerResidualOutput4AD();
+
+    // stop recording
+    this->globalADTape_.setPassive();
+
+    // set seeds to the outputs
+    this->assignSeeds2ResidualGradient(seeds);
+
+    // now calculate the reverse matrix-vector product
+    this->globalADTape_.evaluate();
+
+    // get the matrix-vector product from the inputs
+    for (label i = 0; i < nParameters; i++)
+    {
+        product[i] = parametersArray[i].getGradient();
+    }
+
+    // clean up AD
+    this->globalADTape_.clearAdjoints();
+    this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    for (label i = 0; i < nParameters; i++)
+    {
+        this->globalADTape_.deactivateValue(parametersArray[i]);
+    }
+    for (label i = 0; i < nParameters; i++)
+    {
+        this->setRegressionParameter(i, parametersArray[i]);
+    }
+    this->regressionModelCompute();
+    this->updateStateBoundaryConditions();
+    this->calcResiduals();
+
+    delete[] parametersArray;
+    delete[] volCoordsArray;
+    delete[] statesArray;
+
+#endif
+}
+
 void DASolver::calcdRdXvTPsiAD(
     const Vec xvVec,
     const Vec wVec,
diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H
index 56d5f8c4..4abdc546 100644
--- a/src/adjoint/DASolver/DASolver.H
+++ b/src/adjoint/DASolver/DASolver.H
@@ -33,6 +33,7 @@
 #include "DAField.H"
 #include "DAPartDeriv.H"
 #include "DALinearEqn.H"
+#include "DARegression.H"
 #include "volPointInterpolation.H"
 #include "IOMRFZoneListDF.H"
 #include "interpolateSplineXY.H"
@@ -99,6 +100,9 @@ protected:
     /// DAStateInfo pointer
     autoPtr<DAStateInfo> daStateInfoPtr_;
 
+    /// DARegression pointer
+    autoPtr<DARegression> daRegressionPtr_;
+
     /// the stateInfo_ list from DAStateInfo object
     HashTable<wordList> stateInfo_;
 
@@ -1127,6 +1131,50 @@ public:
         return forwardADDerivVal_[objFuncName];
     }
 
+    /// get the number of regression model parameters
+    label getNRegressionParameters()
+    {
+#ifdef SolidDASolver
+        return 0;
+#else
+        return daRegressionPtr_->nParameters();
+#endif
+    }
+
+    /// get the regression parameter
+    scalar getRegressionParameter(const label idxI)
+    {
+#ifdef SolidDASolver
+        return 0.0;
+#else
+        return daRegressionPtr_->getParameter(idxI);
+#endif
+    }
+
+    /// set the regression parameter
+    void setRegressionParameter(const label idxI, scalar val)
+    {
+#ifndef SolidDASolver
+        daRegressionPtr_->setParameter(idxI, val);
+#endif
+    }
+
+    /// call the compute method of the regression model
+    void regressionModelCompute()
+    {
+#ifndef SolidDASolver
+        daRegressionPtr_->compute();
+#endif
+    }
+
+    /// calculte [dR/dRegParameters]^T * psi
+    void calcdRdRegParTPsiAD(
+        const double* volCoords,
+        const double* states,
+        const double* parameters,
+        const double* seeds,
+        double* product);
+
     /// update the primal state boundary condition based on the primalBC dict
     void setPrimalBoundaryConditions(const label printInfo = 1);
 
diff --git a/src/adjoint/Make/files_Compressible b/src/adjoint/Make/files_Compressible
index 1901dcca..bdd91278 100755
--- a/src/adjoint/Make/files_Compressible
+++ b/src/adjoint/Make/files_Compressible
@@ -26,6 +26,8 @@ DAModel/DARadiationModel/DAP1.C
 
 DAModel/DAModel.C
 
+DARegression/DARegression.C
+
 DAStateInfo/DAStateInfo.C
 DAStateInfo/DAStateInfoRhoSimpleFoam.C
 DAStateInfo/DAStateInfoRhoSimpleCFoam.C
diff --git a/src/adjoint/Make/files_Incompressible b/src/adjoint/Make/files_Incompressible
index aec2dcf2..072f4eb8 100755
--- a/src/adjoint/Make/files_Incompressible
+++ b/src/adjoint/Make/files_Incompressible
@@ -27,6 +27,8 @@ DAModel/DARadiationModel/DAP1.C
 
 DAModel/DAModel.C
 
+DARegression/DARegression.C
+
 DAStateInfo/DAStateInfo.C
 DAStateInfo/DAStateInfoSimpleFoam.C
 DAStateInfo/DAStateInfoPisoFoam.C
diff --git a/src/include/createAdjointCompressible.H b/src/include/createAdjointCompressible.H
index f42a4267..ad7fad91 100755
--- a/src/include/createAdjointCompressible.H
+++ b/src/include/createAdjointCompressible.H
@@ -27,6 +27,8 @@ daTurbulenceModelPtr_.reset(DATurbulenceModel::New(turbModelName, mesh, daOption
 
 daModelPtr_.reset(new DAModel(mesh, daOptionPtr_()));
 
+daRegressionPtr_.reset(new DARegression(mesh, daOptionPtr_(), daModelPtr_()));
+
 word solverName = daOptionPtr_->getOption<word>("solverName");
 daStateInfoPtr_.reset(DAStateInfo::New(solverName, mesh, daOptionPtr_(), daModelPtr_()));
 stateInfo_ = daStateInfoPtr_->getStateInfo();
diff --git a/src/include/createAdjointIncompressible.H b/src/include/createAdjointIncompressible.H
index 065ffb6e..7cc6987c 100755
--- a/src/include/createAdjointIncompressible.H
+++ b/src/include/createAdjointIncompressible.H
@@ -27,6 +27,8 @@ daTurbulenceModelPtr_.reset(DATurbulenceModel::New(turbModelName, mesh, daOption
 
 daModelPtr_.reset(new DAModel(mesh, daOptionPtr_()));
 
+daRegressionPtr_.reset(new DARegression(mesh, daOptionPtr_(), daModelPtr_()));
+
 word solverName = daOptionPtr_->getOption<word>("solverName");
 daStateInfoPtr_.reset(DAStateInfo::New(solverName, mesh, daOptionPtr_(), daModelPtr_()));
 stateInfo_ = daStateInfoPtr_->getStateInfo();
diff --git a/src/include/setForwardADSeeds.H b/src/include/setForwardADSeeds.H
index fd993a0e..7228ad4e 100755
--- a/src/include/setForwardADSeeds.H
+++ b/src/include/setForwardADSeeds.H
@@ -371,6 +371,19 @@ if (daOptionPtr_->getAllOptions().subDict("useAD").getWord("mode") == "forward")
             }
         }
     }
+    else if (type == "RegPar")
+    {
+        // for ACTD derivs
+        scalar seed = 0.0;
+
+        label seedIndex = daOptionPtr_->getAllOptions().subDict("useAD").getLabel("seedIndex");
+
+        seed = daRegressionPtr_->getParameter(seedIndex);
+
+        seed.setGradient(1.0);
+
+        daRegressionPtr_->setParameter(seedIndex, seed);
+    }
 }
 
 #endif
diff --git a/src/pyDASolvers/DASolvers.H b/src/pyDASolvers/DASolvers.H
index 55694ff6..bf280fce 100755
--- a/src/pyDASolvers/DASolvers.H
+++ b/src/pyDASolvers/DASolvers.H
@@ -582,6 +582,36 @@ public:
                                    << abort(FatalError);
     }
 
+    /// get the number of regression model parameters
+    label getNRegressionParameters()
+    {
+        return DASolverPtr_->getNRegressionParameters();
+    }
+
+    /// call the compute method of the regression model
+    void regressionModelCompute()
+    {
+        DASolverPtr_->regressionModelCompute();
+    }
+
+    /// set the regression parameter
+    void setRegressionParameter(label idxI, double val)
+    {
+        scalar val1 = val;
+        DASolverPtr_->setRegressionParameter(idxI, val1);
+    }
+
+    /// calculte [dR/dRegParameters]^T * psi
+    void calcdRdRegParTPsiAD(
+        const double* volCoords,
+        const double* states,
+        const double* parameters,
+        const double* seeds,
+        double* product)
+    {
+        DASolverPtr_->calcdRdRegParTPsiAD(volCoords, states, parameters, seeds, product);
+    }
+
     void calcdRdThermalTPsiAD(
         const double* volCoords,
         const double* states,
diff --git a/src/pyDASolvers/pyDASolvers.pyx b/src/pyDASolvers/pyDASolvers.pyx
index 79916961..055e1933 100755
--- a/src/pyDASolvers/pyDASolvers.pyx
+++ b/src/pyDASolvers/pyDASolvers.pyx
@@ -75,6 +75,7 @@ cdef extern from "DASolvers.H" namespace "Foam":
         void calcdRdFieldTPsiAD(PetscVec, PetscVec, PetscVec, char *, PetscVec)
         void calcdFdFieldAD(PetscVec, PetscVec, char *, char *, PetscVec)
         void calcdRdThermalTPsiAD(double *, double *, double *, double *, double *)
+        void calcdRdRegParTPsiAD(double *, double *, double *, double *, double *)
         void calcdRdWOldTPsiAD(int, PetscVec, PetscVec)
         void convertMPIVec2SeqVec(PetscVec, PetscVec)
         void syncDAOptionToActuatorDVs()
@@ -103,6 +104,9 @@ cdef extern from "DASolvers.H" namespace "Foam":
         void getThermal(double *, double *, double *)
         void getThermalAD(char *, double *, double *, double *, double *)
         void setThermal(double *)
+        int getNRegressionParameters()
+        void setRegressionParameter(int, double)
+        void regressionModelCompute()
         void getOFField(char *, char *, PetscVec)
         void getAcousticData(PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, char*)
         void printAllOptions()
@@ -471,6 +475,41 @@ cdef class pyDASolvers:
         assert len(thermal) == self.getNCouplingFaces() * 2, "invalid array size!"
         cdef double *thermal_data = <double*>thermal.data
         self._thisptr.setThermal(thermal_data)
+    
+    def getNRegressionParameters(self):
+        return self._thisptr.getNRegressionParameters()
+    
+    def setRegressionParameter(self, idx, val):
+        self._thisptr.setRegressionParameter(idx, val)
+    
+    def regressionModelCompute(self):
+        self._thisptr.regressionModelCompute()
+    
+    def calcdRdRegParTPsiAD(self, 
+            np.ndarray[double, ndim=1, mode="c"] volCoords,
+            np.ndarray[double, ndim=1, mode="c"] states,
+            np.ndarray[double, ndim=1, mode="c"] parameters,
+            np.ndarray[double, ndim=1, mode="c"] seeds,
+            np.ndarray[double, ndim=1, mode="c"] product):
+        
+        assert len(volCoords) == self.getNLocalPoints() * 3, "invalid array size!"
+        assert len(states) == self.getNLocalAdjointStates(), "invalid array size!"
+        assert len(parameters) == self.getNRegressionParameters(), "invalid array size!"
+        assert len(seeds) == self.getNLocalAdjointStates(), "invalid array size!"
+        assert len(product) == self.getNRegressionParameters(), "invalid array size!"
+
+        cdef double *volCoords_data = <double*>volCoords.data
+        cdef double *states_data = <double*>states.data
+        cdef double *parameters_data = <double*>parameters.data
+        cdef double *seeds_data = <double*>seeds.data
+        cdef double *product_data = <double*>product.data
+
+        self._thisptr.calcdRdRegParTPsiAD(
+            volCoords_data, 
+            states_data, 
+            parameters_data, 
+            seeds_data, 
+            product_data)
 
     def getAcousticData(self, Vec x, Vec y, Vec z, Vec nX, Vec nY, Vec nZ, Vec a, Vec fX, Vec fY, Vec fZ, groupName):
         self._thisptr.getAcousticData(x.vec, y.vec, z.vec, nX.vec, nY.vec, nZ.vec, a.vec, fX.vec, fY.vec, fZ.vec, groupName)
diff --git a/tests/refs/DAFoam_Test_DAPimpleFoamRef.txt b/tests/refs/DAFoam_Test_DAPimpleFoamRef.txt
index c07c907e..9f5b6a15 100755
--- a/tests/refs/DAFoam_Test_DAPimpleFoamRef.txt
+++ b/tests/refs/DAFoam_Test_DAPimpleFoamRef.txt
@@ -1,27 +1,29 @@
 Dictionary Key: CD
-@value    7.766963355160421 1e-08 1e-10
+@value    7.730954054805606 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
 Dictionary Key: CD
 Dictionary Key: actuator
-@value  0.06183736794552468 0.0001 1e-06
-@value -0.02964620159378894 0.0001 1e-06
-@value 8.590768616224275e-07 0.0001 1e-06
+@value  0.06101115696083378 0.0001 1e-06
+@value -0.02920189888090502 0.0001 1e-06
+@value 8.496815758355939e-07 0.0001 1e-06
 @value                    0 0.0001 1e-06
-@value -0.02866393949934062 0.0001 1e-06
-@value -5.444357559822659e-06 0.0001 1e-06
-@value  -0.1105254610216093 0.0001 1e-06
-@value -0.004255864874662887 0.0001 1e-06
-@value 0.0002947623588585657 0.0001 1e-06
-@value -6.411774379853559e-07 0.0001 1e-06
-@value -0.02604158402942627 0.0001 1e-06
-@value -0.01969385064592553 0.0001 1e-06
+@value -0.02812087805700929 0.0001 1e-06
+@value -5.376812853366829e-06 0.0001 1e-06
+@value  -0.1084638902133479 0.0001 1e-06
+@value -0.003563872996223476 0.0001 1e-06
+@value 0.0002903835775881846 0.0001 1e-06
+@value -5.804106441305712e-07 0.0001 1e-06
+@value -0.02560009058308888 0.0001 1e-06
+@value -0.01945709084472788 0.0001 1e-06
 @value                    0 0.0001 1e-06
 Dictionary Key: alpha
-@value    0.297696232336795 0.0001 1e-06
+@value   0.2967470635402201 0.0001 1e-06
+Dictionary Key: parameter
+@value   0.1569591562230965 0.0001 1e-06
 Dictionary Key: twist
-@value   0.5121846706565929 0.0001 1e-06
+@value    0.507329805800139 0.0001 1e-06
 Dictionary Key: uin
-@value    0.189403455340968 0.0001 1e-06
+@value   0.1879122164504213 0.0001 1e-06
 Dictionary Key: fail
 @value                    0 0.0001 1e-06
\ No newline at end of file
diff --git a/tests/runTests_DAPimpleFoam.py b/tests/runTests_DAPimpleFoam.py
index 46bddcb9..9d344a47 100755
--- a/tests/runTests_DAPimpleFoam.py
+++ b/tests/runTests_DAPimpleFoam.py
@@ -34,7 +34,12 @@
     "designSurfaces": ["wing"],
     "solverName": "DAPimpleFoam",
     "primalBC": {"U0": {"variable": "U", "patches": ["inout"], "value": [U0, 0, 0]}, "useWallFunction": True},
-    "unsteadyAdjoint": {"mode": "timeAccurate", "PCMatPrecomputeInterval": 5, "PCMatUpdateInterval": 1, "objFuncTimeOperator": "average"},
+    "unsteadyAdjoint": {
+        "mode": "timeAccurate",
+        "PCMatPrecomputeInterval": 5,
+        "PCMatUpdateInterval": 1,
+        "objFuncTimeOperator": "average",
+    },
     "printIntervalUnsteady": 1,
     "fvSource": {
         "disk1": {
@@ -54,6 +59,14 @@
             "targetThrust": 1.0,
         },
     },
+    "regressionModel": {
+        "active": True,
+        "modelType": "neuralNetwork",
+        "inputNames": ["PoD", "SoQ", "chiSA", "pGradStream"],
+        "outputName": "betaFI",
+        "hiddenLayerNeurons": [5, 5],
+        "outputShift": 1.0,
+    },
     "objFunc": {
         "CD": {
             "part1": {
@@ -95,7 +108,12 @@
         "uin": {"designVarType": "BC", "patches": ["inout"], "variable": "U", "comp": 0},
         "twist": {"designVarType": "FFD"},
         "alpha": {"designVarType": "AOA", "patches": ["inout"], "flowAxis": "x", "normalAxis": "y"},
-        "actuator": {"actuatorName": "disk1", "designVarType": "ACTD", "comps": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]},
+        "actuator": {
+            "actuatorName": "disk1",
+            "designVarType": "ACTD",
+            "comps": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],
+        },
+        "parameter": {"designVarType": "RegPar"},
     },
 }
 
@@ -188,10 +206,23 @@ def actuator(val, geo):
     scale=1.0,
 )
 
+
+def regModel(val, geo):
+    for idxI in range(len(val)):
+        val1 = float(val[idxI])
+        DASolver.setRegressionParameter(idxI, val1)
+
+
 # =============================================================================
 # DAFoam initialization
 # =============================================================================
 DASolver = PYDAFOAM(options=daOptions, comm=gcomm)
+
+nParameters = DASolver.solver.getNRegressionParameters()
+
+parameter0 = np.ones(nParameters) * 0.1
+DVGeo.addGlobalDV("parameter", parameter0, regModel, lower=-100.0, upper=100.0, scale=1.0)
+
 DASolver.setDVGeo(DVGeo)
 mesh = USMesh(options=meshOptions, comm=gcomm)
 DASolver.printFamilyList()
@@ -226,6 +257,9 @@ def actuator(val, geo):
     funcsSens = {}
     funcsSens, fail = optFuncs.calcObjFuncSens(xDV, funcs)
 
+    parameterNormU = np.linalg.norm(funcsSens["CD"]["parameter"])
+    funcsSens["CD"]["parameter"] = parameterNormU
+
     if gcomm.rank == 0:
         reg_write_dict(funcs, 1e-8, 1e-10)
         reg_write_dict(funcsSens, 1e-4, 1e-6)