From be232cb4e2a7b9c03f32c19f005e749d17e2bb95 Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Sun, 14 Apr 2024 12:44:14 -0500
Subject: [PATCH] Enabled using multiple regModels (#629)

* Enabled using multiple regModels.

* Fixed a warning for location obj.

* Fixed a regModel bug for forward.
---
 dafoam/mphys/mphys_dafoam.py                  |  22 +-
 dafoam/pyDAFoam.py                            |  80 ++-
 src/adjoint/DAObjFunc/DAObjFuncLocation.C     |   2 +-
 src/adjoint/DARegression/DARegression.C       | 481 ++++++++++--------
 src/adjoint/DARegression/DARegression.H       | 111 ++--
 .../DASolver/DAPimpleFoam/DAPimpleFoam.C      |   3 -
 .../DARhoPimpleFoam/DARhoPimpleFoam.C         |   3 -
 .../DARhoSimpleCFoam/DARhoSimpleCFoam.C       |   3 -
 .../DARhoSimpleFoam/DARhoSimpleFoam.C         |   3 -
 .../DASolver/DASimpleFoam/DASimpleFoam.C      |   3 -
 .../DASolver/DASimpleTFoam/DASimpleTFoam.C    |   3 -
 src/adjoint/DASolver/DASolver.C               |  22 +-
 src/adjoint/DASolver/DASolver.H               |  14 +-
 .../DASolver/DATurboFoam/DATurboFoam.C        |   3 -
 src/include/setForwardADSeeds.H               |   6 +-
 src/pyDASolvers/DASolvers.H                   |  14 +-
 src/pyDASolvers/pyDASolvers.pyx               |  28 +-
 .../DAFoam_Test_DARhoSimpleFoamFIMLRef.txt    |  14 +-
 tests/runTests_DAPimpleFoam.py                |  37 +-
 tests/runTests_DARhoSimpleFoamFIML.py         | 156 +++---
 20 files changed, 558 insertions(+), 450 deletions(-)

diff --git a/dafoam/mphys/mphys_dafoam.py b/dafoam/mphys/mphys_dafoam.py
index 417d07f4..51ce0be2 100644
--- a/dafoam/mphys/mphys_dafoam.py
+++ b/dafoam/mphys/mphys_dafoam.py
@@ -664,7 +664,8 @@ def setup(self):
                     distributed = designVariables[dvName]["distributed"]
                 self.add_input(dvName, distributed=distributed, shape_by_conn=True, tags=["mphys_coupling"])
             elif dvType == "RegPar":
-                nParameters = self.DASolver.solver.getNRegressionParameters()
+                modelName = designVariables[dvName]["modelName"]
+                nParameters = self.DASolver.solver.getNRegressionParameters(modelName.encode())
                 self.add_input(dvName, distributed=False, shape=nParameters, tags=["mphys_coupling"])
             else:
                 raise AnalysisError("designVarType %s not supported! " % dvType)
@@ -915,7 +916,11 @@ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
                         product = np.zeros_like(parameters)
                         seeds = d_residuals["%s_states" % self.discipline]
 
-                        DASolver.solverAD.calcdRdRegParTPsiAD(volCoords, states, parameters, seeds, product)
+                        modelName = designVariables[inputName]["modelName"]
+
+                        DASolver.solverAD.calcdRdRegParTPsiAD(
+                            volCoords, states, parameters, seeds, modelName.encode(), product
+                        )
 
                         # reduce to make sure all procs get the same product
                         product = self.comm.allreduce(product, op=MPI.SUM)
@@ -1238,7 +1243,8 @@ def setup(self):
                     distributed = designVariables[dvName]["distributed"]
                 self.add_input(dvName, distributed=distributed, shape_by_conn=True, tags=["mphys_coupling"])
             elif dvType == "RegPar":
-                nParameters = self.DASolver.solver.getNRegressionParameters()
+                modelName = designVariables[dvName]["modelName"]
+                nParameters = self.DASolver.solver.getNRegressionParameters(modelName.encode())
                 self.add_input(dvName, distributed=False, shape=nParameters, tags=["mphys_coupling"])
             else:
                 raise AnalysisError("designVarType %s not supported! " % dvType)
@@ -1456,8 +1462,16 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
                         parameters = inputs[inputName]
                         product = np.zeros_like(parameters)
 
+                        modelName = designVariables[inputName]["modelName"]
+
                         DASolver.solverAD.calcdFdRegParAD(
-                            volCoords, states, parameters, objFuncName.encode(), inputName.encode(), product
+                            volCoords,
+                            states,
+                            parameters,
+                            objFuncName.encode(),
+                            inputName.encode(),
+                            modelName.encode(),
+                            product,
                         )
 
                         d_inputs[inputName] += product * fBar
diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py
index 8be65864..52570295 100755
--- a/dafoam/pyDAFoam.py
+++ b/dafoam/pyDAFoam.py
@@ -514,21 +514,43 @@ def __init__(self):
         self.useConstrainHbyA = False
 
         ## parameters for regression models
+        ## we support defining multiple regression models. Each regression model can have only one output
+        ## but it can have multiple input features. Refer to src/adjoint/DARegression/DARegression.C for
+        ## a full list of supported input features. There are two supported regression model types:
+        ## neural network and radial basis function. We can shift and scale the inputs and outputs
+        ## we can also prescribe a default output value. The default output will be used in resetting
+        ## when they are nan, inf, or out of the prescribe upper and lower bounds
         self.regressionModel = {
             "active": False,
-            "modelType": "neuralNetwork",
-            "inputNames": ["None"],
-            "outputName": "None",
-            "hiddenLayerNeurons": [0],
-            "inputShift": [0.0],
-            "inputScale": [1.0],
-            "outputShift": 0.0,
-            "outputScale": 1.0,
-            "outputUpperBound": 1e8,
-            "outputLowerBound": -1e8,
-            "activationFunction": "sigmoid",
-            "printInputInfo": True,
-            "defaultOutputValue": 0.0,
+            # "model1": {
+            #    "modelType": "neuralNetwork",
+            #    "inputNames": ["PoD", "VoS", "PSoSS", "chiSA"],
+            #    "outputName": "betaFINuTilda",
+            #    "hiddenLayerNeurons": [20, 20],
+            #    "inputShift": [0.0],
+            #    "inputScale": [1.0],
+            #    "outputShift": 1.0,
+            #    "outputScale": 1.0,
+            #    "outputUpperBound": 1e1,
+            #    "outputLowerBound": -1e1,
+            #    "activationFunction": "sigmoid",  # other options are relu and tanh
+            #    "printInputInfo": True,
+            #    "defaultOutputValue": 1.0,
+            # },
+            # "model2": {
+            #    "modelType": "radialBasisFunction",
+            #    "inputNames": ["KoU2", "ReWall", "CoP", "TauoK"],
+            #    "outputName": "betaFIOmega",
+            #    "nRBFs": 50,
+            #    "inputShift": [0.0],
+            #    "inputScale": [1.0],
+            #    "outputShift": 1.0,
+            #    "outputScale": 1.0,
+            #    "outputUpperBound": 1e1,
+            #    "outputLowerBound": -1e1,
+            #    "printInputInfo": True,
+            #    "defaultOutputValue": 1.0,
+            # },
         }
 
         ## whether to use averaged states
@@ -539,10 +561,7 @@ def __init__(self):
         ## we use no time steps for averaging. 0.8 means we use the last 20% of the time step for averaging.
         ## We usually don't use 0 because the flow will need some spin up time, so using the spin-up
         ## flow field for meanStates will be slightly inaccurate.
-        self.useMeanStates = {
-            "active": False,
-            "start": 0.5
-        }
+        self.useMeanStates = {"active": False, "start": 0.5}
 
         # *********************************************************************************************
         # ************************************ Advance Options ****************************************
@@ -2181,7 +2200,7 @@ def calcTotalDerivsACT(self, objFuncName, designVarName, designVarType, dFScalin
             else:
                 self.adjTotalDeriv[objFuncName][designVarName][i] = totalDerivSeq[i]
 
-    def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accumulateTotal=False):
+    def calcTotalDerivsRegPar(self, objFuncName, designVarName, modelName, dFScaling=1.0, accumulateTotal=False):
 
         nDVs = 0
         parameters = None
@@ -2195,7 +2214,7 @@ def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accum
             nDVs = len(self.internalDV[designVarName]["init"])
             parameters = self.internalDV[designVarName]["value"]
 
-        nParameters = self.solver.getNRegressionParameters()
+        nParameters = self.getNRegressionParameters(modelName)
         if nDVs != nParameters:
             raise Error("number of parameters not valid!")
 
@@ -2207,12 +2226,12 @@ def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accum
 
         # calc dFdRegPar
         self.solverAD.calcdFdRegParAD(
-            xvArray, wArray, parameters, objFuncName.encode(), designVarName.encode(), dFdRegPar
+            xvArray, wArray, parameters, objFuncName.encode(), designVarName.encode(), modelName.encode(), dFdRegPar
         )
         dFdRegPar *= dFScaling
 
         # calculate dRdFieldT*Psi and save it to totalDeriv
-        self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, totalDerivArray)
+        self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, modelName.encode(), totalDerivArray)
         # all reduce because parameters is a global DV
         totalDerivArray = self.comm.allreduce(totalDerivArray, op=MPI.SUM)
 
@@ -2437,7 +2456,8 @@ def solveAdjointUnsteady(self):
                             fieldType = designVarDict[designVarName]["fieldType"]
                             self.calcTotalDerivsField(objFuncName, designVarName, fieldType, dFScaling, True)
                         elif designVarDict[designVarName]["designVarType"] == "RegPar":
-                            self.calcTotalDerivsRegPar(objFuncName, designVarName, dFScaling, True)
+                            modelName = designVarDict[designVarName]["modelName"]
+                            self.calcTotalDerivsRegPar(objFuncName, designVarName, modelName, dFScaling, True)
                         else:
                             raise Error("designVarType not valid!")
 
@@ -2736,7 +2756,8 @@ def solveAdjoint(self):
                     # loop over all objectives
                     for objFuncName in objFuncDict:
                         if objFuncName in self.objFuncNames4Adj:
-                            self.calcTotalDerivsRegPar(objFuncName, designVarName)
+                            modelName = designVarDict[designVarName]["modelName"]
+                            self.calcTotalDerivsRegPar(objFuncName, designVarName, modelName)
                 else:
                     raise Error("For Field design variable type, we only support useAD->mode=reverse")
             else:
@@ -4013,7 +4034,7 @@ def setThermal(self, thermalArray):
         if self.getOption("useAD")["mode"] in ["forward", "reverse"]:
             self.solverAD.setThermal(thermalArray)
 
-    def setRegressionParameter(self, idx, val):
+    def setRegressionParameter(self, modelName, idx, val):
         """
         Update the regression parameters
 
@@ -4026,9 +4047,16 @@ def setRegressionParameter(self, idx, val):
             the parameter value to set
         """
 
-        self.solver.setRegressionParameter(idx, val)
+        self.solver.setRegressionParameter(modelName.encode(), idx, val)
         if self.getOption("useAD")["mode"] in ["forward", "reverse"]:
-            self.solverAD.setRegressionParameter(idx, val)
+            self.solverAD.setRegressionParameter(modelName.encode(), idx, val)
+
+    def getNRegressionParameters(self, modelName):
+        """
+        Get the number of regression model parameters
+        """
+        nParameters = self.solver.getNRegressionParameters(modelName.encode())
+        return nParameters
 
     def setFieldValue4GlobalCellI(self, fieldName, val, globalCellI, compI=0):
         """
diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.C b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
index d0fe804e..a568f107 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncLocation.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
@@ -72,7 +72,7 @@ DAObjFuncLocation::DAObjFuncLocation(
         // we assume the patchI and faceI do not change during the optimization
         // otherwise, we should use maxRadiusKS instead
         scalar maxR = -100000;
-        label maxRPatchI, maxRFaceI;
+        label maxRPatchI = -999, maxRFaceI = -999;
         forAll(objFuncFaceSources_, idxI)
         {
             const label& objFuncFaceI = objFuncFaceSources_[idxI];
diff --git a/src/adjoint/DARegression/DARegression.C b/src/adjoint/DARegression/DARegression.C
index c6be2943..5ed93a2c 100755
--- a/src/adjoint/DARegression/DARegression.C
+++ b/src/adjoint/DARegression/DARegression.C
@@ -23,95 +23,121 @@ DARegression::DARegression(
       daModel_(daModel)
 {
     dictionary regSubDict = daOption.getAllOptions().subDict("regressionModel");
-    regSubDict.readEntry<word>("modelType", modelType_);
-    regSubDict.readEntry<wordList>("inputNames", inputNames_);
-    regSubDict.readEntry<word>("outputName", outputName_);
+    active_ = regSubDict.getLabel("active");
 
-    regSubDict.readEntry<scalarList>("inputShift", inputShift_);
+    // initialize parameters
+    if (active_)
+    {
+        forAll(regSubDict.toc(), idxI)
+        {
+            word key = regSubDict.toc()[idxI];
+            if (key != "active")
+            {
+                modelNames_.append(key);
+            }
+        }
 
-    regSubDict.readEntry<scalarList>("inputScale", inputScale_);
+        forAll(modelNames_, idxI)
+        {
+            word modelName = modelNames_[idxI];
+            dictionary modelSubDict = daOption.getAllOptions().subDict("regressionModel").subDict(modelName);
 
-    regSubDict.readEntry<scalar>("outputShift", outputShift_);
+            modelType_.set(modelName, modelSubDict.getWord("modelType"));
 
-    regSubDict.readEntry<scalar>("outputScale", outputScale_);
+            wordList tempWordList;
+            scalarList tempScalarList;
+            labelList tempLabelList;
 
-    regSubDict.readEntry<scalar>("outputUpperBound", outputUpperBound_);
+            modelSubDict.readEntry<wordList>("inputNames", tempWordList);
+            inputNames_.set(modelName, tempWordList);
 
-    regSubDict.readEntry<scalar>("outputLowerBound", outputLowerBound_);
+            outputName_.set(modelName, modelSubDict.getWord("outputName"));
 
-    regSubDict.readEntry<label>("printInputInfo", printInputInfo_);
+            modelSubDict.readEntry<scalarList>("inputShift", tempScalarList);
+            inputShift_.set(modelName, tempScalarList);
 
-    regSubDict.readEntry<scalar>("defaultOutputValue", defaultOutputValue_);
+            modelSubDict.readEntry<scalarList>("inputScale", tempScalarList);
+            inputScale_.set(modelName, tempScalarList);
 
-    writeFeatures_ = regSubDict.lookupOrDefault<label>("writeFeatures", 0);
+            outputShift_.set(modelName, modelSubDict.getScalar("outputShift"));
 
-    active_ = regSubDict.getLabel("active");
+            outputScale_.set(modelName, modelSubDict.getScalar("outputScale"));
 
-    // initialize parameters
-    if (active_)
-    {
-        if (modelType_ == "neuralNetwork")
-        {
-            regSubDict.readEntry<labelList>("hiddenLayerNeurons", hiddenLayerNeurons_);
-            regSubDict.readEntry<word>("activationFunction", activationFunction_);
-            if (activationFunction_ == "ReLU")
+            outputUpperBound_.set(modelName, modelSubDict.getScalar("outputUpperBound"));
+
+            outputLowerBound_.set(modelName, modelSubDict.getScalar("outputLowerBound"));
+
+            defaultOutputValue_.set(modelName, modelSubDict.getScalar("defaultOutputValue"));
+
+            printInputInfo_.set(modelName, modelSubDict.getLabel("printInputInfo"));
+
+            writeFeatures_.set(modelName, modelSubDict.lookupOrDefault<label>("writeFeatures", 0));
+
+            if (modelType_[modelName] == "neuralNetwork")
             {
-                leakyCoeff_ = regSubDict.lookupOrDefault<scalar>("leakyCoeff", 0.0);
+                modelSubDict.readEntry<labelList>("hiddenLayerNeurons", tempLabelList);
+                hiddenLayerNeurons_.set(modelName, tempLabelList);
+                activationFunction_.set(modelName, modelSubDict.getWord("activationFunction"));
+                if (activationFunction_[modelName] == "relu")
+                {
+                    leakyCoeff_.set(modelName, modelSubDict.lookupOrDefault<scalar>("leakyCoeff", 0.0));
+                }
+            }
+            else if (modelType_[modelName] == "radialBasisFunction")
+            {
+                nRBFs_.set(modelName, modelSubDict.getLabel("nRBFs"));
+            }
+            else
+            {
+                FatalErrorIn("DARegression") << "modelType_: " << modelType_[modelName] << " not supported. Options are: neuralNetwork and radialBasisFunction" << abort(FatalError);
             }
-        }
-        else if (modelType_ == "radialBasisFunction")
-        {
-            nRBFs_ = regSubDict.getLabel("nRBFs");
-        }
-        else
-        {
-            FatalErrorIn("DARegression") << "modelType_: " << modelType_ << " not supported. Options are: neuralNetwork and radialBasisFunction" << abort(FatalError);
-        }
 
-        // check the sizes
-        if (inputNames_.size() != inputShift_.size() || inputNames_.size() != inputScale_.size())
-        {
-            FatalErrorIn("DARegression") << "inputNames has different sizes than inputShift or inputScale" << abort(FatalError);
-        }
+            // check the sizes
+            if (inputNames_[modelName].size() != inputShift_[modelName].size()
+                || inputNames_[modelName].size() != inputScale_[modelName].size())
+            {
+                FatalErrorIn("DARegression") << "inputNames has different sizes than inputShift or inputScale" << abort(FatalError);
+            }
 
-        label nParameters = this->nParameters();
-        parameters_.setSize(nParameters);
-        forAll(parameters_, idxI)
-        {
-            parameters_[idxI] = 0.0;
-        }
+            label nParameters = this->nParameters(modelName);
+            scalarList parameters(nParameters, 0.0);
+            parameters_.set(modelName, parameters);
 
-        // initialize the ptr scalarFields
-        features_.setSize(inputNames_.size());
-        forAll(features_, idxI)
-        {
-            word inputName = inputNames_[idxI];
-
-            features_.set(
-                idxI,
-                new volScalarField(
-                    IOobject(
-                        inputName,
-                        mesh_.time().timeName(),
+            // initialize the ptr scalarFields
+            label nInputs = inputNames_[modelName].size();
+            PtrList<volScalarField> features(nInputs);
+            forAll(inputNames_[modelName], idxI)
+            {
+                word inputName = inputNames_[modelName][idxI];
+                features.set(
+                    idxI,
+                    new volScalarField(
+                        IOobject(
+                            inputName,
+                            mesh_.time().timeName(),
+                            mesh_,
+                            IOobject::NO_READ,
+                            IOobject::NO_WRITE),
                         mesh_,
-                        IOobject::NO_READ,
-                        IOobject::NO_WRITE),
-                    mesh_,
-                    dimensionedScalar(inputName, dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0)));
+                        dimensionedScalar(inputName, dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0)));
+            }
+            features_.set(modelName, features);
         }
     }
 }
 
-void DARegression::calcInputFeatures()
+void DARegression::calcInputFeatures(word modelName)
 {
     /*
     Description:
-        Calculate the input features
+        Calculate the input features. Here features is a unique list for all the regModel's inputNames.
+        This is because two regModel's inputNames can have overlap, so we don't want to compute
+        duplicated features
     */
 
-    forAll(inputNames_, idxI)
+    forAll(features_[modelName], idxI)
     {
-        word inputName = inputNames_[idxI];
+        word inputName = inputNames_[modelName][idxI];
         if (inputName == "VoS")
         {
             // vorticity / strain
@@ -120,21 +146,21 @@ void DARegression::calcInputFeatures()
             const volTensorField& gradU = tgradU();
             volScalarField magOmega = mag(skew(gradU));
             volScalarField magS = mag(symm(gradU));
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
-                features_[idxI][cellI] = (magOmega[cellI] / (magS[cellI] + 1e-16) + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (magOmega[cellI] / (magS[cellI] + 1e-16) + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else if (inputName == "PoD")
         {
             // production / destruction
-            daModel_.getTurbProdOverDestruct(features_[idxI]);
-            forAll(features_[idxI], cellI)
+            daModel_.getTurbProdOverDestruct(features_[modelName][idxI]);
+            forAll(features_[modelName][idxI], cellI)
             {
-                features_[idxI][cellI] = (features_[idxI][cellI] + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (features_[modelName][idxI][cellI] + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else if (inputName == "chiSA")
         {
@@ -142,11 +168,11 @@ void DARegression::calcInputFeatures()
             // the chi() function from SA
             const volScalarField& nuTilda = mesh_.thisDb().lookupObject<volScalarField>("nuTilda");
             volScalarField nu = daModel_.getDATurbulenceModel().nu();
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
-                features_[idxI][cellI] = (nuTilda[cellI] / nu[cellI] + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (nuTilda[cellI] / nu[cellI] + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
 #endif
         }
         else if (inputName == "pGradStream")
@@ -162,11 +188,11 @@ void DARegression::calcInputFeatures()
             }
             volScalarField::Internal pGradAlongStream = (U() & pGrad()) / pG_denominator();
 
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
-                features_[idxI][cellI] = (pGradAlongStream[cellI] + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (pGradAlongStream[cellI] + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else if (inputName == "PSoSS")
         {
@@ -178,15 +204,15 @@ void DARegression::calcInputFeatures()
             volVectorField pGrad("gradP", fvc::grad(p));
             vector diagUGrad = vector::zero;
             scalar val = 0;
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
                 diagUGrad[0] = gradU[cellI].xx();
                 diagUGrad[1] = gradU[cellI].yy();
                 diagUGrad[2] = gradU[cellI].zz();
                 val = mag(pGrad[cellI]) / (mag(pGrad[cellI]) + mag(3.0 * cmptAv(U[cellI] & diagUGrad)) + 1e-16);
-                features_[idxI][cellI] = (val + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (val + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else if (inputName == "SCurv")
         {
@@ -196,12 +222,12 @@ void DARegression::calcInputFeatures()
             const volTensorField& gradU = tgradU();
 
             scalar val = 0;
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
                 val = mag(U[cellI] & gradU[cellI]) / (mag(U[cellI] & U[cellI]) + mag(U[cellI] & gradU[cellI]) + 1e-16);
-                features_[idxI][cellI] = (val + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (val + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else if (inputName == "UOrth")
         {
@@ -211,12 +237,12 @@ void DARegression::calcInputFeatures()
             const volTensorField& gradU = tgradU();
 
             scalar val = 0;
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
                 val = mag(U[cellI] & gradU[cellI] & U[cellI]) / (mag(U[cellI]) * mag(gradU[cellI] & U[cellI]) + mag(U[cellI] & gradU[cellI] & U[cellI]) + 1e-16);
-                features_[idxI][cellI] = (val + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (val + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else if (inputName == "KoU2")
         {
@@ -224,12 +250,12 @@ void DARegression::calcInputFeatures()
             const volScalarField& k = mesh_.thisDb().lookupObject<volScalarField>("k");
             const volVectorField& U = mesh_.thisDb().lookupObject<volVectorField>("U");
             scalar val = 0;
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
                 val = k[cellI] / (0.5 * (U[cellI] & U[cellI]) + 1e-16);
-                features_[idxI][cellI] = (val + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (val + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else if (inputName == "ReWall")
         {
@@ -239,23 +265,23 @@ void DARegression::calcInputFeatures()
             const volScalarField& k = mesh_.thisDb().lookupObject<volScalarField>("k");
             volScalarField nu = daModel_.getDATurbulenceModel().nu();
             scalar val = 0;
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
                 val = sqrt(k[cellI]) * y[cellI] / (50.0 * nu[cellI]);
-                features_[idxI][cellI] = (val + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (val + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
 #endif
         }
         else if (inputName == "CoP")
         {
             // convective / production
-            daModel_.getTurbConvOverProd(features_[idxI]);
-            forAll(features_[idxI], cellI)
+            daModel_.getTurbConvOverProd(features_[modelName][idxI]);
+            forAll(features_[modelName][idxI], cellI)
             {
-                features_[idxI][cellI] = (features_[idxI][cellI] + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (features_[modelName][idxI][cellI] + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else if (inputName == "TauoK")
         {
@@ -265,12 +291,12 @@ void DARegression::calcInputFeatures()
             const volVectorField& U = mesh_.thisDb().lookupObject<volVectorField>("U");
             volSymmTensorField tau(2.0 / 3.0 * I * k - nut * twoSymm(fvc::grad(U)));
             scalar val = 0;
-            forAll(features_[idxI], cellI)
+            forAll(features_[modelName][idxI], cellI)
             {
                 val = mag(tau[cellI]) / (k[cellI] + 1e-16);
-                features_[idxI][cellI] = (val + inputShift_[idxI]) * inputScale_[idxI];
+                features_[modelName][idxI][cellI] = (val + inputShift_[modelName][idxI]) * inputScale_[modelName][idxI];
             }
-            features_[idxI].correctBoundaryConditions();
+            features_[modelName][idxI].correctBoundaryConditions();
         }
         else
         {
@@ -304,151 +330,158 @@ label DARegression::compute()
         a volScalarField prescribed by outputName
     */
 
-    // if the output variable is not found in the Db, just return and do nothing
-    if (!active_ || !mesh_.thisDb().foundObject<volScalarField>(outputName_))
+    if (!active_)
     {
         return 0;
     }
 
     label fail = 0;
 
-    volScalarField& outputField = const_cast<volScalarField&>(mesh_.thisDb().lookupObject<volScalarField>(outputName_));
-
-    if (modelType_ == "neuralNetwork")
+    forAll(modelNames_, idxI)
     {
-        label nHiddenLayers = hiddenLayerNeurons_.size();
-
-        this->calcInputFeatures();
+        word modelName = modelNames_[idxI];
 
-        List<List<scalar>> layerVals;
-        layerVals.setSize(nHiddenLayers);
-        for (label layerI = 0; layerI < nHiddenLayers; layerI++)
+        // if the output variable is not found in the Db, just return and do nothing
+        if (!mesh_.thisDb().foundObject<volScalarField>(outputName_[modelName]))
         {
-            label nNeurons = hiddenLayerNeurons_[layerI];
-            layerVals[layerI].setSize(nNeurons);
+            return 0;
         }
 
-        forAll(mesh_.cells(), cellI)
-        {
-            label counterI = 0;
+        // compute the inputFeature for all inputs
+        this->calcInputFeatures(modelName);
+
+        volScalarField& outputField = const_cast<volScalarField&>(mesh_.thisDb().lookupObject<volScalarField>(outputName_[modelName]));
 
+        if (modelType_[modelName] == "neuralNetwork")
+        {
+            label nHiddenLayers = hiddenLayerNeurons_[modelName].size();
+            List<List<scalar>> layerVals;
+            layerVals.setSize(nHiddenLayers);
             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++)
+                label nNeurons = hiddenLayerNeurons_[modelName][layerI];
+                layerVals[layerI].setSize(nNeurons);
+            }
+
+            forAll(mesh_.cells(), cellI)
+            {
+                label counterI = 0;
+
+                for (label layerI = 0; layerI < nHiddenLayers; layerI++)
                 {
-                    if (layerI == 0)
+                    label nNeurons = hiddenLayerNeurons_[modelName][layerI];
+                    forAll(layerVals[layerI], neuronI)
                     {
-                        // for the 1st hidden layer, we use the input layer as the input
-                        forAll(inputNames_, neuronJ)
-                        {
-                            // weighted sum
-                            layerVals[layerI][neuronI] += features_[neuronJ][cellI] * parameters_[counterI];
-                            counterI++;
-                        }
+                        layerVals[layerI][neuronI] = 0.0;
                     }
-                    else
+                    for (label neuronI = 0; neuronI < nNeurons; neuronI++)
                     {
-                        // for the rest of hidden layer, we use the previous hidden layer as the input
-                        forAll(layerVals[layerI - 1], neuronJ)
+                        if (layerI == 0)
                         {
-                            // weighted sum
-                            layerVals[layerI][neuronI] += layerVals[layerI - 1][neuronJ] * parameters_[counterI];
-                            counterI++;
+                            // for the 1st hidden layer, we use the input layer as the input
+                            forAll(inputNames_[modelName], neuronJ)
+                            {
+                                // weighted sum
+                                layerVals[layerI][neuronI] += features_[modelName][neuronJ][cellI] * parameters_[modelName][counterI];
+                                counterI++;
+                            }
                         }
-                    }
-                    // bias
-                    layerVals[layerI][neuronI] += parameters_[counterI];
-                    counterI++;
-                    // activation function
-                    if (activationFunction_ == "sigmoid")
-                    {
-                        layerVals[layerI][neuronI] = 1 / (1 + exp(-layerVals[layerI][neuronI]));
-                    }
-                    else if (activationFunction_ == "tanh")
-                    {
-                        layerVals[layerI][neuronI] = (1 - exp(-2 * layerVals[layerI][neuronI])) / (1 + exp(-2 * layerVals[layerI][neuronI]));
-                    }
-                    else if (activationFunction_ == "ReLU")
-                    {
-                        if (layerVals[layerI][neuronI] < 0)
+                        else
                         {
-                            layerVals[layerI][neuronI] = leakyCoeff_ * layerVals[layerI][neuronI];
+                            // for the rest of hidden layer, we use the previous hidden layer as the input
+                            forAll(layerVals[layerI - 1], neuronJ)
+                            {
+                                // weighted sum
+                                layerVals[layerI][neuronI] += layerVals[layerI - 1][neuronJ] * parameters_[modelName][counterI];
+                                counterI++;
+                            }
+                        }
+                        // bias
+                        layerVals[layerI][neuronI] += parameters_[modelName][counterI];
+                        counterI++;
+                        // activation function
+                        if (activationFunction_[modelName] == "sigmoid")
+                        {
+                            layerVals[layerI][neuronI] = 1 / (1 + exp(-layerVals[layerI][neuronI]));
+                        }
+                        else if (activationFunction_[modelName] == "tanh")
+                        {
+                            layerVals[layerI][neuronI] = (1 - exp(-2 * layerVals[layerI][neuronI])) / (1 + exp(-2 * layerVals[layerI][neuronI]));
+                        }
+                        else if (activationFunction_[modelName] == "relu")
+                        {
+                            if (layerVals[layerI][neuronI] < 0)
+                            {
+                                layerVals[layerI][neuronI] = leakyCoeff_[modelName] * layerVals[layerI][neuronI];
+                            }
+                        }
+                        else
+                        {
+                            FatalErrorIn("") << "activationFunction not valid. Options are: sigmoid, tanh, and relu" << abort(FatalError);
                         }
-                    }
-                    else
-                    {
-                        FatalErrorIn("") << "activationFunction not valid. Options are: sigmoid, tanh, and ReLU" << abort(FatalError);
                     }
                 }
-            }
-            // 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];
-
-            // no activation function for the output layer
-
-            outputField[cellI] = outputScale_ * (outputVal + outputShift_);
-        }
+                // 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_[modelName][counterI];
+                    counterI++;
+                }
+                // bias
+                outputVal += parameters_[modelName][counterI];
 
-        // check if the output values are valid otherwise fix/bound them
-        fail = this->checkOutput(outputField);
+                // no activation function for the output layer
 
-        outputField.correctBoundaryConditions();
-    }
-    else if (modelType_ == "radialBasisFunction")
-    {
+                outputField[cellI] = outputScale_[modelName] * (outputVal + outputShift_[modelName]);
+            }
 
-        this->calcInputFeatures();
+            // check if the output values are valid otherwise fix/bound them
+            fail += this->checkOutput(modelName, outputField);
 
-        label nInputs = inputNames_.size();
+            outputField.correctBoundaryConditions();
+        }
+        else if (modelType_[modelName] == "radialBasisFunction")
+        {
+            label nInputs = inputNames_[modelName].size();
 
-        // increment of the parameters for each RBF basis
-        label dP = 2 * nInputs + 1;
+            // increment of the parameters for each RBF basis
+            label dP = 2 * nInputs + 1;
 
-        forAll(mesh_.cells(), cellI)
-        {
-            scalar outputVal = 0.0;
-            for (label i = 0; i < nRBFs_; i++)
+            forAll(mesh_.cells(), cellI)
             {
-                scalar expCoeff = 0.0;
-                for (label j = 0; j < nInputs; j++)
+                scalar outputVal = 0.0;
+                for (label i = 0; i < nRBFs_[modelName]; i++)
                 {
-                    scalar A = (features_[j][cellI] - parameters_[dP * i + 2 * j]) * (features_[j][cellI] - parameters_[dP * i + 2 * j]);
-                    scalar B = 2 * parameters_[dP * i + 2 * j + 1] * parameters_[dP * i + 2 * j + 1];
-                    expCoeff += A / B;
+                    scalar expCoeff = 0.0;
+                    for (label j = 0; j < nInputs; j++)
+                    {
+                        scalar A = (features_[modelName][j][cellI] - parameters_[modelName][dP * i + 2 * j]) * (features_[modelName][j][cellI] - parameters_[modelName][dP * i + 2 * j]);
+                        scalar B = 2 * parameters_[modelName][dP * i + 2 * j + 1] * parameters_[modelName][dP * i + 2 * j + 1];
+                        expCoeff += A / B;
+                    }
+                    outputVal += parameters_[modelName][dP * i + dP - 1] * exp(-expCoeff);
                 }
-                outputVal += parameters_[dP * i + dP - 1] * exp(-expCoeff);
-            }
 
-            outputField[cellI] = outputScale_ * (outputVal + outputShift_);
-        }
+                outputField[cellI] = outputScale_[modelName] * (outputVal + outputShift_[modelName]);
+            }
 
-        // check if the output values are valid otherwise fix/bound them
-        fail = this->checkOutput(outputField);
+            // check if the output values are valid otherwise fix/bound them
+            fail += this->checkOutput(modelName, outputField);
 
-        outputField.correctBoundaryConditions();
-    }
-    else
-    {
-        FatalErrorIn("") << "modelType_: " << modelType_ << " not supported. Options are: neuralNetwork and radialBasisFunction" << abort(FatalError);
+            outputField.correctBoundaryConditions();
+        }
+        else
+        {
+            FatalErrorIn("") << "modelType_: " << modelType_ << " not supported. Options are: neuralNetwork and radialBasisFunction" << abort(FatalError);
+        }
     }
 
     return fail;
 }
 
-label DARegression::nParameters()
+label DARegression::nParameters(word modelName)
 {
     /*
     Description:
@@ -460,49 +493,49 @@ label DARegression::nParameters()
         FatalErrorIn("") << "nParameters() is called but the regression model is not active!" << abort(FatalError);
     }
 
-    if (modelType_ == "neuralNetwork")
+    if (modelType_[modelName] == "neuralNetwork")
     {
-        label nHiddenLayers = hiddenLayerNeurons_.size();
-        label nInputs = inputNames_.size();
+        label nHiddenLayers = hiddenLayerNeurons_[modelName].size();
+        label nInputs = inputNames_[modelName].size();
 
         // add weights
         // input
-        label nParameters = nInputs * hiddenLayerNeurons_[0];
+        label nParameters = nInputs * hiddenLayerNeurons_[modelName][0];
         // hidden layers
-        for (label layerI = 1; layerI < hiddenLayerNeurons_.size(); layerI++)
+        for (label layerI = 1; layerI < hiddenLayerNeurons_[modelName].size(); layerI++)
         {
-            nParameters += hiddenLayerNeurons_[layerI] * hiddenLayerNeurons_[layerI - 1];
+            nParameters += hiddenLayerNeurons_[modelName][layerI] * hiddenLayerNeurons_[modelName][layerI - 1];
         }
         // output
-        nParameters += hiddenLayerNeurons_[nHiddenLayers - 1] * 1;
+        nParameters += hiddenLayerNeurons_[modelName][nHiddenLayers - 1] * 1;
 
         // add biases
         // add hidden layers
-        for (label layerI = 0; layerI < hiddenLayerNeurons_.size(); layerI++)
+        for (label layerI = 0; layerI < hiddenLayerNeurons_[modelName].size(); layerI++)
         {
-            nParameters += hiddenLayerNeurons_[layerI];
+            nParameters += hiddenLayerNeurons_[modelName][layerI];
         }
         // add output layer
         nParameters += 1;
 
         return nParameters;
     }
-    else if (modelType_ == "radialBasisFunction")
+    else if (modelType_[modelName] == "radialBasisFunction")
     {
-        label nInputs = inputNames_.size();
+        label nInputs = inputNames_[modelName].size();
 
         // each RBF has a weight, nInputs mean, and nInputs std
-        label nParameters = nRBFs_ * (2 * nInputs + 1);
+        label nParameters = nRBFs_[modelName] * (2 * nInputs + 1);
 
         return nParameters;
     }
     else
     {
-        FatalErrorIn("") << "modelType_: " << modelType_ << " not supported. Options are: neuralNetwork and radialBasisFunction" << abort(FatalError);
+        FatalErrorIn("") << "modelType_: " << modelType_[modelName] << " not supported. Options are: neuralNetwork and radialBasisFunction" << abort(FatalError);
     }
 }
 
-label DARegression::checkOutput(volScalarField& outputField)
+label DARegression::checkOutput(word modelName, volScalarField& outputField)
 {
     /*
     Description:
@@ -523,40 +556,40 @@ label DARegression::checkOutput(volScalarField& outputField)
     {
         if (std::isnan(outputField[cellI]))
         {
-            outputField[cellI] = defaultOutputValue_;
+            outputField[cellI] = defaultOutputValue_[modelName];
             isNaN = 1;
         }
         if (std::isinf(outputField[cellI]))
         {
-            outputField[cellI] = defaultOutputValue_;
+            outputField[cellI] = defaultOutputValue_[modelName];
             isInf = 1;
         }
-        if (outputField[cellI] > outputUpperBound_)
+        if (outputField[cellI] > outputUpperBound_[modelName])
         {
-            outputField[cellI] = outputUpperBound_;
+            outputField[cellI] = outputUpperBound_[modelName];
             isBounded = 1;
         }
-        if (outputField[cellI] < outputLowerBound_)
+        if (outputField[cellI] < outputLowerBound_[modelName])
         {
-            outputField[cellI] = outputLowerBound_;
+            outputField[cellI] = outputLowerBound_[modelName];
             isBounded = 1;
         }
     }
     if (isBounded == 1)
     {
-        Pout << "************* Warning! output values are bounded between " << outputLowerBound_ << " and " << outputUpperBound_ << endl;
+        Pout << "************* Warning! output values are bounded between " << outputLowerBound_[modelName] << " and " << outputUpperBound_[modelName] << endl;
         fail = 1;
     }
 
     if (isNaN == 1)
     {
-        Pout << "************* Warning! output values have nan and are set to " << defaultOutputValue_ << endl;
+        Pout << "************* Warning! output values have nan and are set to " << defaultOutputValue_[modelName] << endl;
         fail = 1;
     }
 
     if (isInf == 1)
     {
-        Pout << "************* Warning! output values have inf and are set to " << defaultOutputValue_ << endl;
+        Pout << "************* Warning! output values have inf and are set to " << defaultOutputValue_[modelName] << endl;
         fail = 1;
     }
 
diff --git a/src/adjoint/DARegression/DARegression.H b/src/adjoint/DARegression/DARegression.H
index 91292f5a..7144daa0 100755
--- a/src/adjoint/DARegression/DARegression.H
+++ b/src/adjoint/DARegression/DARegression.H
@@ -48,62 +48,65 @@ protected:
     /// DAModel object
     const DAModel& daModel_;
 
+    /// whether the regression model is active
+    label active_;
+
+    /// name of the regression models
+    wordList modelNames_;
+
     /// the type of regression model
-    word modelType_;
+    HashTable<word> modelType_;
 
     /// a list of words for the inputs
-    wordList inputNames_;
+    HashTable<wordList> inputNames_;
 
     /// a list of words for the outputs
-    word outputName_;
+    HashTable<word> outputName_;
 
     /// number of neurons hidden layers of the neural network
-    labelList hiddenLayerNeurons_;
+    HashTable<labelList> hiddenLayerNeurons_;
 
     /// we can shift each input. we always shift before scaling it.
-    scalarList inputShift_;
+    HashTable<scalarList> inputShift_;
 
     /// we can scale each input. we always shift before scaling it.
-    scalarList inputScale_;
+    HashTable<scalarList> inputScale_;
 
     /// we can shift the output. we always shift before scaling it.
-    scalar outputShift_;
+    HashTable<scalar> outputShift_;
 
     /// we can scale the output. we always shift before scaling it.
-    scalar outputScale_;
+    HashTable<scalar> outputScale_;
 
     /// the parameters for the regression model
-    scalarList parameters_;
-
-    /// whether the regression model is active
-    label active_;
+    HashTable<scalarList> parameters_;
 
     /// neural network activation function
-    word activationFunction_;
+    HashTable<word> activationFunction_;
 
     /// if the ReLU activation function is used we can prescribe a potentially leaky coefficient
-    scalar leakyCoeff_ = 0.0;
+    HashTable<scalar> leakyCoeff_;
 
     /// the upper bound for the output
-    scalar outputUpperBound_;
+    HashTable<scalar> outputUpperBound_;
 
     /// the lower bound for the output
-    scalar outputLowerBound_;
+    HashTable<scalar> outputLowerBound_;
 
     /// whether to print the input range info this is used to scale the input
-    label printInputInfo_;
+    HashTable<label> printInputInfo_;
 
     /// default output values
-    scalar defaultOutputValue_;
+    HashTable<scalar> defaultOutputValue_;
 
     /// number of radial basis function
-    label nRBFs_;
+    HashTable<label> nRBFs_;
 
     /// a list of scalarField to save the input features
-    PtrList<volScalarField> features_;
+    HashTable<PtrList<volScalarField>> features_;
 
     /// whether to write the feature fields to the disk
-    label writeFeatures_ = 0;
+    HashTable<label> writeFeatures_;
 
 public:
     /// Constructors
@@ -123,55 +126,45 @@ public:
     label compute();
 
     /// calculate the input flow features
-    void calcInputFeatures();
+    void calcInputFeatures(word modelName);
 
     /// get the number of parameters for this regression model
-    label nParameters();
+    label nParameters(word modelName);
 
     /// get a specific parameter value
-    scalar getParameter(label idxI)
+    scalar getParameter(word modelName, label idxI)
     {
-        return parameters_[idxI];
+        return parameters_[modelName][idxI];
     }
 
     /// set a value to a parameter give an index and a value
-    void setParameter(label idxI, scalar val)
+    void setParameter(word modelName, 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);
-            }
-        }
+        parameters_[modelName][idxI] = val;
     }
 
     /// check if the output values are valid otherwise bound or fix them
-    label checkOutput(volScalarField& outputField);
+    label checkOutput(word modelName, volScalarField& outputField);
 
     /// print the input
     void printInputInfo()
     {
-        if (active_ && printInputInfo_)
+        if (!active_)
         {
-            this->calcInputFeatures();
-
-            forAll(inputNames_, idxI)
+            return;
+        }
+        forAll(modelNames_, idxI)
+        {
+            word modelName = modelNames_[idxI];
+            if (printInputInfo_[modelName])
             {
-                word name = inputNames_[idxI];
-                Info << "RegModel input info. " << name << " Min: " << gMin(features_[idxI]) << " Max: " << gMax(features_[idxI])
-                     << " Avg: " << gAverage(features_[idxI]) << " Std: " << sqrt(gSumSqr(features_[idxI]) / features_[idxI].size()) << endl;
+                Info << "RegModel input info for " << modelName << endl;
+                forAll(inputNames_[modelName], idxI)
+                {
+                    word name = inputNames_[modelName][idxI];
+                    Info << name << " Min: " << gMin(features_[modelName][idxI]) << " Max: " << gMax(features_[modelName][idxI])
+                         << " Avg: " << gAverage(features_[modelName][idxI]) << " Std: " << sqrt(gSumSqr(features_[modelName][idxI]) / features_[modelName][idxI].size()) << endl;
+                }
             }
         }
     }
@@ -179,11 +172,19 @@ public:
     /// write the features to the disk
     void writeFeatures()
     {
-        if (writeFeatures_)
+        if (!active_)
+        {
+            return;
+        }
+        forAll(modelNames_, idxI)
         {
-            forAll(features_, idxI)
+            word modelName = modelNames_[idxI];
+            if (writeFeatures_[modelName])
             {
-                features_[idxI].write();
+                forAll(features_[modelName], idxI)
+                {
+                    features_[modelName][idxI].write();
+                }
             }
         }
     }
diff --git a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C
index 1a479842..1186bf99 100755
--- a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C
+++ b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C
@@ -173,9 +173,6 @@ 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
     label regModelFail = 0;
     label fail = 0;
diff --git a/src/adjoint/DASolver/DARhoPimpleFoam/DARhoPimpleFoam.C b/src/adjoint/DASolver/DARhoPimpleFoam/DARhoPimpleFoam.C
index 624f0c81..b53f1fad 100755
--- a/src/adjoint/DASolver/DARhoPimpleFoam/DARhoPimpleFoam.C
+++ b/src/adjoint/DASolver/DARhoPimpleFoam/DARhoPimpleFoam.C
@@ -172,9 +172,6 @@ label DARhoPimpleFoam::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
     label regModelFail = 0;
     label fail = 0;
diff --git a/src/adjoint/DASolver/DARhoSimpleCFoam/DARhoSimpleCFoam.C b/src/adjoint/DASolver/DARhoSimpleCFoam/DARhoSimpleCFoam.C
index 08702961..22bdd6cf 100755
--- a/src/adjoint/DASolver/DARhoSimpleCFoam/DARhoSimpleCFoam.C
+++ b/src/adjoint/DASolver/DARhoSimpleCFoam/DARhoSimpleCFoam.C
@@ -129,9 +129,6 @@ label DARhoSimpleCFoam::solvePrimal(
         }
     }
 
-    // check if the parameters are set in the Python layer
-    daRegressionPtr_->validate();
-
     // if useMeanStates is used, we need to zero meanStates before the primal run
     this->zeroMeanStates();
     
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
index 1d492985..3627ebac 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
@@ -146,9 +146,6 @@ label DARhoSimpleFoam::solvePrimal(
         }
     }
 
-    // check if the parameters are set in the Python layer
-    daRegressionPtr_->validate();
-
     // if useMeanStates is used, we need to zero meanStates before the primal run
     this->zeroMeanStates();
 
diff --git a/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C b/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
index 000c46ed..30c72ee7 100755
--- a/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
+++ b/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
@@ -156,9 +156,6 @@ label DASimpleFoam::solvePrimal(
         }
     }
 
-    // check if the parameters are set in the Python layer
-    daRegressionPtr_->validate();
-
     // if useMeanStates is used, we need to zero meanStates before the primal run
     this->zeroMeanStates();
 
diff --git a/src/adjoint/DASolver/DASimpleTFoam/DASimpleTFoam.C b/src/adjoint/DASolver/DASimpleTFoam/DASimpleTFoam.C
index f40a1b53..57fd0c3b 100755
--- a/src/adjoint/DASolver/DASimpleTFoam/DASimpleTFoam.C
+++ b/src/adjoint/DASolver/DASimpleTFoam/DASimpleTFoam.C
@@ -136,9 +136,6 @@ label DASimpleTFoam::solvePrimal(
     // if the forwardModeAD is active, we need to set the seed here
 #include "setForwardADSeeds.H"
 
-    // check if the parameters are set in the Python layer
-    daRegressionPtr_->validate();
-
     // if useMeanStates is used, we need to zero meanStates before the primal run
     this->zeroMeanStates();
 
diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
index 9bc608f7..0289b1ed 100644
--- a/src/adjoint/DASolver/DASolver.C
+++ b/src/adjoint/DASolver/DASolver.C
@@ -5900,6 +5900,7 @@ void DASolver::calcdFdRegParAD(
     const double* parameters,
     const word objFuncName,
     const word designVarName,
+    const word modelName,
     double* dFdRegPar)
 {
 #ifdef CODI_AD_REVERSE
@@ -5916,12 +5917,14 @@ void DASolver::calcdFdRegParAD(
         objFuncName: the name of the objective function
 
         designVarName: name of the design variable
+
+        modelName: regression model name defined in the regressionModel dict
     
     Output:
         dFdRegPar: dF/dRegPar
     */
 
-    Info << "Calculating dFdRegPar using reverse-mode AD" << endl;
+    Info << "Calculating dFdRegPar using reverse-mode AD: " << modelName << endl;
 
     scalar* volCoordsArray = new scalar[daIndexPtr_->nLocalXv];
     for (label i = 0; i < daIndexPtr_->nLocalXv; i++)
@@ -5935,7 +5938,7 @@ void DASolver::calcdFdRegParAD(
         statesArray[i] = states[i];
     }
 
-    label nParameters = this->getNRegressionParameters();
+    label nParameters = this->getNRegressionParameters(modelName);
     scalar* parametersArray = new scalar[nParameters];
     for (label i = 0; i < nParameters; i++)
     {
@@ -5974,7 +5977,7 @@ void DASolver::calcdFdRegParAD(
         }
         for (label i = 0; i < nParameters; i++)
         {
-            this->setRegressionParameter(i, parametersArray[i]);
+            this->setRegressionParameter(modelName, i, parametersArray[i]);
         }
 
         // update the BC, this func will also call regressionModelCompute
@@ -6026,7 +6029,7 @@ void DASolver::calcdFdRegParAD(
         }
         for (label i = 0; i < nParameters; i++)
         {
-            this->setRegressionParameter(i, parametersArray[i]);
+            this->setRegressionParameter(modelName, i, parametersArray[i]);
         }
 
         this->updateStateBoundaryConditions();
@@ -6160,6 +6163,7 @@ void DASolver::calcdRdRegParTPsiAD(
     const double* states,
     const double* parameters,
     const double* seeds,
+    const word modelName,
     double* product)
 {
 #ifdef CODI_AD_REVERSE
@@ -6176,12 +6180,14 @@ void DASolver::calcdRdRegParTPsiAD(
         parameters: the parameters for the regression model
 
         psi: the vector to multiply dRdXv
+
+        modelName: the regression model name defined in the regressionModel dict
     
     Output:
         prodVec: the matrix-vector products [dR/dRegParameters]^T*Psi
     */
 
-    Info << "Calculating [dRdRegPar]^T * Psi using reverse-mode AD" << endl;
+    Info << "Calculating [dRdRegPar]^T * Psi using reverse-mode AD: " << modelName << endl;
 
     scalar* volCoordsArray = new scalar[daIndexPtr_->nLocalXv];
     for (label i = 0; i < daIndexPtr_->nLocalXv; i++)
@@ -6195,7 +6201,7 @@ void DASolver::calcdRdRegParTPsiAD(
         statesArray[i] = states[i];
     }
 
-    label nParameters = this->getNRegressionParameters();
+    label nParameters = this->getNRegressionParameters(modelName);
     scalar* parametersArray = new scalar[nParameters];
     for (label i = 0; i < nParameters; i++)
     {
@@ -6219,7 +6225,7 @@ void DASolver::calcdRdRegParTPsiAD(
     }
     for (label i = 0; i < nParameters; i++)
     {
-        this->setRegressionParameter(i, parametersArray[i]);
+        this->setRegressionParameter(modelName, i, parametersArray[i]);
     }
 
     // update the BC, this func will also call regressionModelCompute
@@ -6260,7 +6266,7 @@ void DASolver::calcdRdRegParTPsiAD(
     }
     for (label i = 0; i < nParameters; i++)
     {
-        this->setRegressionParameter(i, parametersArray[i]);
+        this->setRegressionParameter(modelName, i, parametersArray[i]);
     }
     this->updateStateBoundaryConditions();
     this->calcResiduals();
diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H
index 7ca2630c..2d69a6b6 100644
--- a/src/adjoint/DASolver/DASolver.H
+++ b/src/adjoint/DASolver/DASolver.H
@@ -612,6 +612,7 @@ public:
         const double* parameters,
         const word objFuncName,
         const word designVarName,
+        const word modelName,
         double* dFdRegPar);
 
     /// compute dRdWOld^T*Psi
@@ -1177,30 +1178,30 @@ public:
     }
 
     /// get the number of regression model parameters
-    label getNRegressionParameters()
+    label getNRegressionParameters(word modelName)
     {
 #ifdef SolidDASolver
         return 0;
 #else
-        return daRegressionPtr_->nParameters();
+        return daRegressionPtr_->nParameters(modelName);
 #endif
     }
 
     /// get the regression parameter
-    scalar getRegressionParameter(const label idxI)
+    scalar getRegressionParameter(word modelName, const label idxI)
     {
 #ifdef SolidDASolver
         return 0.0;
 #else
-        return daRegressionPtr_->getParameter(idxI);
+        return daRegressionPtr_->getParameter(modelName, idxI);
 #endif
     }
 
     /// set the regression parameter
-    void setRegressionParameter(const label idxI, scalar val)
+    void setRegressionParameter(word modelName, const label idxI, scalar val)
     {
 #ifndef SolidDASolver
-        daRegressionPtr_->setParameter(idxI, val);
+        daRegressionPtr_->setParameter(modelName, idxI, val);
 #endif
     }
 
@@ -1218,6 +1219,7 @@ public:
         const double* states,
         const double* parameters,
         const double* seeds,
+        const word modelName,
         double* product);
 
     /// update the primal state boundary condition based on the primalBC dict
diff --git a/src/adjoint/DASolver/DATurboFoam/DATurboFoam.C b/src/adjoint/DASolver/DATurboFoam/DATurboFoam.C
index faacc390..0d247c99 100755
--- a/src/adjoint/DASolver/DATurboFoam/DATurboFoam.C
+++ b/src/adjoint/DASolver/DATurboFoam/DATurboFoam.C
@@ -99,9 +99,6 @@ label DATurboFoam::solvePrimal(
     // if the forwardModeAD is active, we need to set the seed here
 #include "setForwardADSeeds.H"
 
-    // check if the parameters are set in the Python layer
-    daRegressionPtr_->validate();
-
     // if useMeanStates is used, we need to zero meanStates before the primal run
     this->zeroMeanStates();
 
diff --git a/src/include/setForwardADSeeds.H b/src/include/setForwardADSeeds.H
index 7228ad4e..2fc7459f 100755
--- a/src/include/setForwardADSeeds.H
+++ b/src/include/setForwardADSeeds.H
@@ -376,13 +376,15 @@ if (daOptionPtr_->getAllOptions().subDict("useAD").getWord("mode") == "forward")
         // for ACTD derivs
         scalar seed = 0.0;
 
+        word modelName = dvSubDict.getWord("modelName");
+
         label seedIndex = daOptionPtr_->getAllOptions().subDict("useAD").getLabel("seedIndex");
 
-        seed = daRegressionPtr_->getParameter(seedIndex);
+        seed = daRegressionPtr_->getParameter(modelName, seedIndex);
 
         seed.setGradient(1.0);
 
-        daRegressionPtr_->setParameter(seedIndex, seed);
+        daRegressionPtr_->setParameter(modelName, seedIndex, seed);
     }
 }
 
diff --git a/src/pyDASolvers/DASolvers.H b/src/pyDASolvers/DASolvers.H
index 02124fc1..689e6c9c 100755
--- a/src/pyDASolvers/DASolvers.H
+++ b/src/pyDASolvers/DASolvers.H
@@ -586,9 +586,9 @@ public:
     }
 
     /// get the number of regression model parameters
-    label getNRegressionParameters()
+    label getNRegressionParameters(word modelName)
     {
-        return DASolverPtr_->getNRegressionParameters();
+        return DASolverPtr_->getNRegressionParameters(modelName);
     }
 
     /// call the compute method of the regression model
@@ -598,10 +598,10 @@ public:
     }
 
     /// set the regression parameter
-    void setRegressionParameter(label idxI, double val)
+    void setRegressionParameter(word modelName, label idxI, double val)
     {
         scalar val1 = val;
-        DASolverPtr_->setRegressionParameter(idxI, val1);
+        DASolverPtr_->setRegressionParameter(modelName, idxI, val1);
     }
 
     /// calculte [dR/dRegParameters]^T * psi
@@ -610,9 +610,10 @@ public:
         const double* states,
         const double* parameters,
         const double* seeds,
+        const word modelName,
         double* product)
     {
-        DASolverPtr_->calcdRdRegParTPsiAD(volCoords, states, parameters, seeds, product);
+        DASolverPtr_->calcdRdRegParTPsiAD(volCoords, states, parameters, seeds, modelName, product);
     }
 
     /// compute dFdRegPar AD
@@ -622,9 +623,10 @@ public:
         const double* parameters,
         const word objFuncName,
         const word designVarName,
+        const word modelName,
         double* dFdRegPar)
     {
-        DASolverPtr_->calcdFdRegParAD(volCoords, states, parameters, objFuncName, designVarName, dFdRegPar);
+        DASolverPtr_->calcdFdRegParAD(volCoords, states, parameters, objFuncName, designVarName, modelName, dFdRegPar);
     }
 
     void calcdRdThermalTPsiAD(
diff --git a/src/pyDASolvers/pyDASolvers.pyx b/src/pyDASolvers/pyDASolvers.pyx
index f74d475b..b030181f 100755
--- a/src/pyDASolvers/pyDASolvers.pyx
+++ b/src/pyDASolvers/pyDASolvers.pyx
@@ -75,8 +75,8 @@ 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 calcdFdRegParAD(double *, double *, double *, char *, char *, double *)
+        void calcdRdRegParTPsiAD(double *, double *, double *, double *, char *, double *)
+        void calcdFdRegParAD(double *, double *, double *, char *, char *, char *, double *)
         void calcdRdWOldTPsiAD(int, PetscVec, PetscVec)
         void convertMPIVec2SeqVec(PetscVec, PetscVec)
         void syncDAOptionToActuatorDVs()
@@ -105,8 +105,8 @@ 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)
+        int getNRegressionParameters(char *)
+        void setRegressionParameter(char *, int, double)
         void regressionModelCompute()
         void getOFField(char *, char *, PetscVec)
         void getAcousticData(PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, char*)
@@ -481,11 +481,11 @@ cdef class pyDASolvers:
         cdef double *thermal_data = <double*>thermal.data
         self._thisptr.setThermal(thermal_data)
     
-    def getNRegressionParameters(self):
-        return self._thisptr.getNRegressionParameters()
+    def getNRegressionParameters(self, modelName):
+        return self._thisptr.getNRegressionParameters(modelName)
     
-    def setRegressionParameter(self, idx, val):
-        self._thisptr.setRegressionParameter(idx, val)
+    def setRegressionParameter(self, modelName, idx, val):
+        self._thisptr.setRegressionParameter(modelName, idx, val)
     
     def regressionModelCompute(self):
         self._thisptr.regressionModelCompute()
@@ -495,13 +495,14 @@ cdef class pyDASolvers:
             np.ndarray[double, ndim=1, mode="c"] states,
             np.ndarray[double, ndim=1, mode="c"] parameters,
             np.ndarray[double, ndim=1, mode="c"] seeds,
+            modelName,
             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(parameters) == self.getNRegressionParameters(modelName), "invalid array size!"
         assert len(seeds) == self.getNLocalAdjointStates(), "invalid array size!"
-        assert len(product) == self.getNRegressionParameters(), "invalid array size!"
+        assert len(product) == self.getNRegressionParameters(modelName), "invalid array size!"
 
         cdef double *volCoords_data = <double*>volCoords.data
         cdef double *states_data = <double*>states.data
@@ -514,6 +515,7 @@ cdef class pyDASolvers:
             states_data, 
             parameters_data, 
             seeds_data, 
+            modelName,
             product_data)
     
     def calcdFdRegParAD(self, 
@@ -522,12 +524,13 @@ cdef class pyDASolvers:
             np.ndarray[double, ndim=1, mode="c"] parameters,
             objFuncName,
             designVarName,
+            modelName,
             np.ndarray[double, ndim=1, mode="c"] dFdRegPar):
         
         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(dFdRegPar) == self.getNRegressionParameters(), "invalid array size!"
+        assert len(parameters) == self.getNRegressionParameters(modelName), "invalid array size!"
+        assert len(dFdRegPar) == self.getNRegressionParameters(modelName), "invalid array size!"
 
         cdef double *volCoords_data = <double*>volCoords.data
         cdef double *states_data = <double*>states.data
@@ -540,6 +543,7 @@ cdef class pyDASolvers:
             parameters_data, 
             objFuncName,
             designVarName,
+            modelName,
             dFdRegPar_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):
diff --git a/tests/refs/DAFoam_Test_DARhoSimpleFoamFIMLRef.txt b/tests/refs/DAFoam_Test_DARhoSimpleFoamFIMLRef.txt
index 85c4abbe..5486f5e1 100755
--- a/tests/refs/DAFoam_Test_DARhoSimpleFoamFIMLRef.txt
+++ b/tests/refs/DAFoam_Test_DARhoSimpleFoamFIMLRef.txt
@@ -1,13 +1,13 @@
 Dictionary Key: VAR
-@value    19.99003057277397 1e-08 1e-10
-Dictionary Key: VAR_RBF
-@value    19.98999618847414 1e-08 1e-10
-Dictionary Key: VAR_ReLU
-@value    19.99092771152761 1e-08 1e-10
+@value    19.98990109204323 1e-08 1e-10
+Dictionary Key: VAR_RBF_ReLU
+@value    19.99022358408163 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
 Dictionary Key: VAR
-Dictionary Key: parameter
-@value 0.002477420843246395 0.0001 1e-06
+Dictionary Key: parameter1
+@value 0.002543197508377707 0.0001 1e-06
+Dictionary Key: parameter2
+@value 0.005268429858563708 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 dbcd9964..754f383d 100755
--- a/tests/runTests_DAPimpleFoam.py
+++ b/tests/runTests_DAPimpleFoam.py
@@ -62,17 +62,22 @@
     },
     "regressionModel": {
         "active": True,
-        "modelType": "neuralNetwork",
-        "inputNames": ["VoS", "PoD", "chiSA", "pGradStream", "PSoSS", "SCurv", "UOrth"],
-        "outputName": "betaFINuTilda",
-        "hiddenLayerNeurons": [5, 5],
-        "inputShift": [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
-        "inputScale": [1.0, 0.00001, 0.01, 1.0, 1.0, 1.0, 1.0],
-        "outputShift": 1.0,
-        "outputScale": 1.0,
-        "activationFunction": "tanh",
-        "printInputInfo": False,
-        "writeFeatures": True,
+        "model1": {
+            "modelType": "neuralNetwork",
+            "inputNames": ["VoS", "PoD", "chiSA", "pGradStream", "PSoSS", "SCurv", "UOrth"],
+            "outputName": "betaFINuTilda",
+            "hiddenLayerNeurons": [5, 5],
+            "inputShift": [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+            "inputScale": [1.0, 0.00001, 0.01, 1.0, 1.0, 1.0, 1.0],
+            "outputShift": 1.0,
+            "outputScale": 1.0,
+            "activationFunction": "tanh",
+            "printInputInfo": True,
+            "writeFeatures": True,
+            "outputUpperBound": 1e2,
+            "outputLowerBound": -1e2,
+            "defaultOutputValue": 1.0,
+        },
     },
     "objFunc": {
         "CD": {
@@ -107,9 +112,9 @@
                 "scale": 1.0,
                 "addToAdjoint": True,
                 "calcRefVar": True,
-                "ref": [0.1, 0.05, 0.04, 0.02, 0.02, 0.01, 0.0, -0.01, -0.01, -0.02]
+                "ref": [0.1, 0.05, 0.04, 0.02, 0.02, 0.01, 0.0, -0.01, -0.01, -0.02],
             }
-        }
+        },
     },
     "adjStateOrdering": "cell",
     "adjEqnOption": {
@@ -133,7 +138,7 @@
             "designVarType": "ACTD",
             "comps": [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],
         },
-        "parameter": {"designVarType": "RegPar"},
+        "parameter": {"designVarType": "RegPar", "modelName": "model1"},
     },
 }
 
@@ -230,7 +235,7 @@ def actuator(val, geo):
 def regModel(val, geo):
     for idxI in range(len(val)):
         val1 = float(val[idxI])
-        DASolver.setRegressionParameter(idxI, val1)
+        DASolver.setRegressionParameter("model1", idxI, val1)
 
 
 # =============================================================================
@@ -238,7 +243,7 @@ def regModel(val, geo):
 # =============================================================================
 DASolver = PYDAFOAM(options=daOptions, comm=gcomm)
 
-nParameters = DASolver.solver.getNRegressionParameters()
+nParameters = DASolver.getNRegressionParameters("model1")
 
 parameter0 = np.ones(nParameters) * 0.1
 DASolver.addInternalDV("parameter", parameter0, regModel, lower=-100.0, upper=100.0, scale=1.0)
diff --git a/tests/runTests_DARhoSimpleFoamFIML.py b/tests/runTests_DARhoSimpleFoamFIML.py
index c937fb21..ae12dd3d 100755
--- a/tests/runTests_DARhoSimpleFoamFIML.py
+++ b/tests/runTests_DARhoSimpleFoamFIML.py
@@ -60,18 +60,38 @@
     },
     "regressionModel": {
         "active": True,
-        "modelType": "neuralNetwork",
-        "inputNames": ["KoU2", "ReWall", "CoP", "TauoK"],
-        "outputName": "betaFIOmega",
-        "hiddenLayerNeurons": [10, 10],
-        "inputShift": [0.0, 0.0, 0.0, 0.0],
-        "inputScale": [1.0, 1.0, 0.0001, 1.0],
-        "outputShift": 1.0,
-        "outputScale": 1.0,
-        "activationFunction": "tanh",
-        "printInputInfo": False,
-        "defaultOutputValue": 1.0,
-        "writeFeatures": True,
+        "model1": {
+            "modelType": "neuralNetwork",
+            "inputNames": ["KoU2", "ReWall", "CoP", "TauoK"],
+            "outputName": "betaFIOmega",
+            "hiddenLayerNeurons": [10, 10],
+            "inputShift": [0.0, 0.0, 0.0, 0.0],
+            "inputScale": [10000.0, 0.1, 0.0001, 1.0],
+            "outputShift": 1.0,
+            "outputScale": 1.0,
+            "activationFunction": "tanh",
+            "printInputInfo": True,
+            "defaultOutputValue": 1.0,
+            "writeFeatures": True,
+            "outputUpperBound": 1e2,
+            "outputLowerBound": -1e2,
+        },
+        "model2": {
+            "modelType": "neuralNetwork",
+            "inputNames": ["KoU2", "PSoSS", "PoD", "TauoK"],
+            "outputName": "betaFIK",
+            "hiddenLayerNeurons": [20, 20],
+            "inputShift": [0.0, 0.0, 0.0, 0.0],
+            "inputScale": [10000.0, 1.0, 0.01, 1.0],
+            "outputShift": 1.0,
+            "outputScale": 1.0,
+            "activationFunction": "tanh",
+            "printInputInfo": True,
+            "defaultOutputValue": 1.0,
+            "writeFeatures": True,
+            "outputUpperBound": 1e2,
+            "outputLowerBound": -1e2,
+        },
     },
     "objFunc": {
         "VAR": {
@@ -96,22 +116,34 @@
     "adjEqnOption": {"gmresRelTol": 1.0e-10, "gmresAbsTol": 1.0e-15, "pcFillLevel": 1, "jacMatReOrdering": "rcm"},
     # Design variable setup
     "designVar": {
-        "parameter": {"designVarType": "RegPar"},
+        "parameter1": {"designVarType": "RegPar", "modelName": "model1"},
+        "parameter2": {"designVarType": "RegPar", "modelName": "model2"},
     },
 }
 
 
-def regModel(val, DASolver):
+def regModel1(val, DASolver):
     for idxI in range(len(val)):
         val1 = float(val[idxI])
-        DASolver.setRegressionParameter(idxI, val1)
+        DASolver.setRegressionParameter("model1", idxI, val1)
+
+
+def regModel2(val, DASolver):
+    for idxI in range(len(val)):
+        val1 = float(val[idxI])
+        DASolver.setRegressionParameter("model2", idxI, val1)
 
 
 # DAFoam
 DASolver = PYDAFOAM(options=aeroOptions, comm=gcomm)
-nParameters = DASolver.solver.getNRegressionParameters()
-parameter0 = np.ones(nParameters) * 0.01
-DASolver.addInternalDV("parameter", parameter0, regModel, lower=-10, upper=10, scale=100.0)
+
+nParameters1 = DASolver.getNRegressionParameters("model1")
+parameter01 = np.ones(nParameters1) * 0.01
+DASolver.addInternalDV("parameter1", parameter01, regModel1, lower=-10, upper=10, scale=100.0)
+
+nParameters2 = DASolver.getNRegressionParameters("model2")
+parameter02 = np.ones(nParameters2) * 0.015
+DASolver.addInternalDV("parameter2", parameter02, regModel2, lower=-10, upper=10, scale=100.0)
 # set evalFuncs
 evalFuncs = []
 DASolver.setEvalFuncs(evalFuncs)
@@ -135,59 +167,59 @@ def regModel(val, DASolver):
     funcsSens = {}
     funcsSens, fail = optFuncs.calcObjFuncSens(iDV, funcs)
 
-    norm = np.linalg.norm(funcsSens["VAR"]["parameter"])
-    funcsSens["VAR"]["parameter"] = norm
+    norm1 = np.linalg.norm(funcsSens["VAR"]["parameter1"])
+    funcsSens["VAR"]["parameter1"] = norm1
+
+    norm2 = np.linalg.norm(funcsSens["VAR"]["parameter2"])
+    funcsSens["VAR"]["parameter2"] = norm2
 
-    # test RBF regression
+    # test RBF regression and ReLU
     aeroOptions["regressionModel"] = {
         "active": True,
-        "modelType": "radialBasisFunction",
-        "inputNames": ["KoU2", "ReWall", "CoP", "TauoK"],
-        "outputName": "betaFIOmega",
-        "nRBFs": 20,
-        "inputShift": [0.0, 0.0, 0.0, 0.0],
-        "inputScale": [1.0, 1.0, 0.0001, 1.0],
-        "outputShift": 1.0,
-        "outputScale": 1.0,
-        "printInputInfo": True,
-        "defaultOutputValue": 1.0,
-        "writeFeatures": True,
+        "model1": {
+            "modelType": "radialBasisFunction",
+            "inputNames": ["KoU2", "ReWall", "CoP", "TauoK"],
+            "outputName": "betaFIOmega",
+            "nRBFs": 20,
+            "inputShift": [0.0, 0.0, 0.0, 0.0],
+            "inputScale": [10000.0, 0.1, 0.0001, 1.0],
+            "outputShift": 1.0,
+            "outputScale": 1.0,
+            "printInputInfo": True,
+            "defaultOutputValue": 1.0,
+            "writeFeatures": True,
+            "outputUpperBound": 1e2,
+            "outputLowerBound": -1e2,
+        },
+        "model2": {
+            "modelType": "neuralNetwork",
+            "inputNames": ["KoU2", "ReWall", "CoP", "TauoK"],
+            "outputName": "betaFIOmega",
+            "hiddenLayerNeurons": [10, 10],
+            "inputShift": [0.0, 0.0, 0.0, 0.0],
+            "inputScale": [10000.0, 0.1, 0.0001, 1.0],
+            "outputShift": 1.0,
+            "outputScale": 1.0,
+            "activationFunction": "relu",
+            "leakyCoeff": 0.1,
+            "printInputInfo": False,
+            "defaultOutputValue": 1.0,
+            "outputUpperBound": 1e2,
+            "outputLowerBound": -1e2,
+        },
     }
 
     DASolver = PYDAFOAM(options=aeroOptions, comm=gcomm)
-    nParameters = DASolver.solver.getNRegressionParameters()
-    parameter0 = np.ones(nParameters) * 0.03
-    DASolver.addInternalDV("parameter", parameter0, regModel, lower=-10, upper=10, scale=100.0)
+    nParameters1 = DASolver.getNRegressionParameters("model1")
+    parameter01 = np.ones(nParameters1) * 0.05
+    DASolver.addInternalDV("parameter1", parameter01, regModel1, lower=-10, upper=10, scale=100.0)
+    nParameters2 = DASolver.getNRegressionParameters("model2")
+    parameter02 = np.ones(nParameters2) * 0.05
+    DASolver.addInternalDV("parameter2", parameter02, regModel2, lower=-10, upper=10, scale=100.0)
     DASolver()
     funcs1 = {}
     DASolver.evalFunctions(funcs1, evalFuncs=["VAR"])
-    funcs["VAR_RBF"] = funcs1["VAR"]
-
-    # test ReLU
-    aeroOptions["regressionModel"] = {
-        "active": True,
-        "modelType": "neuralNetwork",
-        "inputNames": ["KoU2", "ReWall", "CoP", "TauoK"],
-        "outputName": "betaFIOmega",
-        "hiddenLayerNeurons": [10, 10],
-        "inputShift": [0.0, 0.0, 0.0, 0.0],
-        "inputScale": [1.0, 1.0, 0.0001, 1.0],
-        "outputShift": 1.0,
-        "outputScale": 1.0,
-        "activationFunction": "ReLU",
-        "leakyCoeff": 0.1,
-        "printInputInfo": False,
-        "defaultOutputValue": 1.0,
-    }
-
-    DASolver = PYDAFOAM(options=aeroOptions, comm=gcomm)
-    nParameters = DASolver.solver.getNRegressionParameters()
-    parameter0 = np.ones(nParameters) * 0.05
-    DASolver.addInternalDV("parameter", parameter0, regModel, lower=-10, upper=10, scale=100.0)
-    DASolver()
-    funcs2 = {}
-    DASolver.evalFunctions(funcs2, evalFuncs=["VAR"])
-    funcs["VAR_ReLU"] = funcs2["VAR"]
+    funcs["VAR_RBF_ReLU"] = funcs1["VAR"]
 
     if gcomm.rank == 0:
         reg_write_dict(funcs, 1e-8, 1e-10)