diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py
index 0df81e9b..ad57c291 100755
--- a/dafoam/pyDAFoam.py
+++ b/dafoam/pyDAFoam.py
@@ -2962,6 +2962,47 @@ def getAcousticData(self, groupName=None):
         # Finally map the vector as required.
         return positions, normals, areas, forces
 
+    def getOFField(self, fieldName, fieldType, distributed=False):
+        """
+        Get the OpenFOAM field variable in mesh.Db() and return it as a serial or parallel array
+
+        Input:
+        ------
+        fieldName: name of the OF field
+
+        fieldType: type of the OF Field, can be either scalar or vector
+
+        distributed: whether the array is distributed (parallel)
+
+        Output:
+        ------
+        field: numpy array that saves the field
+
+        """
+
+        if fieldType == "scalar":
+            fieldComp = 1
+        elif fieldType == "vector":
+            fieldComp = 3
+        else:
+            raise Error("fieldType not valid!")
+        nLocalCells = self.solver.getNLocalCells()
+
+        vecSize = fieldComp * nLocalCells
+        fieldVec = PETSc.Vec().create(comm=PETSc.COMM_WORLD)
+        fieldVec.setSizes((vecSize, PETSc.DECIDE), bsize=1)
+        fieldVec.setFromOptions()
+        fieldVec.zeroEntries()
+
+        self.solver.getOFField(fieldName.encode(), fieldType.encode(), fieldVec)
+
+        if distributed:
+            field = self.vec2Array(fieldVec)
+        else:
+            field = self.convertMPIVec2SeqArray(fieldVec)
+
+        return field
+
     def calcTotalDeriv(self, dRdX, dFdX, psi, totalDeriv):
         """
         Compute total derivative
diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
index 3917da31..838a2e92 100644
--- a/src/adjoint/DASolver/DASolver.C
+++ b/src/adjoint/DASolver/DASolver.C
@@ -1130,6 +1130,49 @@ void DASolver::getThermalAD(
 #endif
 }
 
+void DASolver::getOFField(
+    const word fieldName,
+    const word fieldType,
+    Vec field) const
+{
+    /*
+    Description:
+        assign a OpenFoam layer field variable in mesh.Db() to field
+    */
+
+    PetscScalar* vecArray;
+    VecGetArray(field, &vecArray);
+
+    if (fieldType == "scalar")
+    {
+        const volScalarField& field = meshPtr_->thisDb().lookupObject<volScalarField>(fieldName);
+        forAll(field, cellI)
+        {
+            assignValueCheckAD(vecArray[cellI], field[cellI]);
+        }
+    }
+    else if (fieldType == "vector")
+    {
+        const volVectorField& field = meshPtr_->thisDb().lookupObject<volVectorField>(fieldName);
+        label localIdx = 0;
+        forAll(field, cellI)
+        {
+            for (label comp = 0; comp < 3; comp++)
+            {
+                assignValueCheckAD(vecArray[localIdx], field[cellI][comp]);
+                localIdx++;
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn("getField") << " fieldType not valid. Options: scalar or vector"
+                                 << abort(FatalError);
+    }
+
+    VecRestoreArray(field, &vecArray);
+}
+
 void DASolver::getForces(Vec fX, Vec fY, Vec fZ)
 {
     /*
diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H
index a6c32ec8..23d7f4d7 100644
--- a/src/adjoint/DASolver/DASolver.H
+++ b/src/adjoint/DASolver/DASolver.H
@@ -689,6 +689,12 @@ public:
         daFieldPtr_->resVec2OFResField(resVec);
     }
 
+    /// get a variable from OF layer
+    void getOFField(
+        const word fieldName,
+        const word fieldType,
+        Vec field) const;
+
     /// write the matrix in binary format
     void writeMatrixBinary(
         const Mat matIn,
diff --git a/src/pyDASolvers/DASolvers.H b/src/pyDASolvers/DASolvers.H
index 9ba311ad..4c5a3c2d 100755
--- a/src/pyDASolvers/DASolvers.H
+++ b/src/pyDASolvers/DASolvers.H
@@ -576,6 +576,15 @@ public:
         DASolverPtr_->calcdRdThermalTPsiAD(volCoords, states, thermal, seeds, product);
     }
 
+    /// get a variable from OF layer
+    void getOFField(
+        const word fieldName,
+        const word fieldType,
+        Vec field)
+    {
+        DASolverPtr_->getOFField(fieldName, fieldType, field);
+    }
+
     /// return the acoustic data values
     void getAcousticData(Vec x, Vec y, Vec z, Vec nX, Vec nY, Vec nZ, Vec a, Vec fX, Vec fY, Vec fZ, word groupName)
     {
diff --git a/src/pyDASolvers/pyDASolvers.pyx b/src/pyDASolvers/pyDASolvers.pyx
index a872b580..9784b472 100755
--- a/src/pyDASolvers/pyDASolvers.pyx
+++ b/src/pyDASolvers/pyDASolvers.pyx
@@ -101,6 +101,7 @@ cdef extern from "DASolvers.H" namespace "Foam":
         void getThermal(double *, double *, double *)
         void getThermalAD(char *, double *, double *, double *, double *)
         void setThermal(double *)
+        void getOFField(char *, char *, PetscVec)
         void getAcousticData(PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, char*)
         void printAllOptions()
         void updateDAOption(object)
@@ -408,6 +409,9 @@ cdef class pyDASolvers:
     def getForces(self, Vec fX, Vec fY, Vec fZ):
         self._thisptr.getForces(fX.vec, fY.vec, fZ.vec)
     
+    def getOFField(self, fieldName, fieldType, Vec fieldVec):
+        self._thisptr.getOFField(fieldName, fieldType, fieldVec.vec)
+    
     def getThermal(self, 
             np.ndarray[double, ndim=1, mode="c"] volCoords,
             np.ndarray[double, ndim=1, mode="c"] states,
diff --git a/tests/runTests_DAPimpleFoamField.py b/tests/runTests_DAPimpleFoamField.py
index 9224d921..08cacc9d 100755
--- a/tests/runTests_DAPimpleFoamField.py
+++ b/tests/runTests_DAPimpleFoamField.py
@@ -28,7 +28,6 @@
 
 twist0 = 30
 U0 = 10
-nCells = 4032
 alpha0 = 0
 
 # test incompressible solvers
@@ -85,17 +84,6 @@
 DVGeo = DVGeometry("./FFD/FFD.xyz")
 DVGeo.addRefAxis("bodyAxis", xFraction=0.25, alignIndex="k")
 
-
-def betaFieldInversion(val, geo):
-    for idxI, v in enumerate(val):
-        DASolver.setFieldValue4GlobalCellI(b"betaFI", v, idxI)
-        DASolver.updateBoundaryConditions(b"betaFI", b"scalar")
-
-
-beta0 = np.ones(nCells, dtype="d")
-beta0[0] = 1.0
-DVGeo.addGlobalDV("beta", value=beta0, func=betaFieldInversion, lower=1e-5, upper=10.0, scale=1.0)
-
 # =============================================================================
 # DAFoam initialization
 # =============================================================================
@@ -107,6 +95,16 @@ def betaFieldInversion(val, geo):
 evalFuncs = []
 DASolver.setEvalFuncs(evalFuncs)
 
+
+def betaFieldInversion(val, geo):
+    for idxI, v in enumerate(val):
+        DASolver.setFieldValue4GlobalCellI(b"betaFI", v, idxI)
+        DASolver.updateBoundaryConditions(b"betaFI", b"scalar")
+
+
+beta0 = DASolver.getOFField("betaFI", "scalar", False)
+DVGeo.addGlobalDV("beta", value=beta0, func=betaFieldInversion, lower=1e-5, upper=10.0, scale=1.0)
+
 # =============================================================================
 # Constraint setup
 # =============================================================================