diff --git a/src/petsc4Foam/solvers/petscSolver.C b/src/petsc4Foam/solvers/petscSolver.C
index 7eafe4a37d757fa45e902ff16dcf6e35a1ce6e57..8b4f4d7a9e9f04487cf148b24d2121b9d1a4cef8 100644
--- a/src/petsc4Foam/solvers/petscSolver.C
+++ b/src/petsc4Foam/solvers/petscSolver.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
     Copyright (C) 2019-2020 Simone Bna
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ License
 #include "fvMesh.H"
 #include "fvMatrices.H"
 #include "globalIndex.H"
+#include "PrecisionAdaptor.H"
 #include "cyclicLduInterface.H"
 #include "cyclicAMILduInterface.H"
 #include "addToRunTimeSelectionTable.H"
@@ -105,10 +106,10 @@ Foam::petscSolver::petscSolver
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::solverPerformance Foam::petscSolver::solve
+Foam::solverPerformance Foam::petscSolver::scalarSolve
 (
-    scalarField& psi,
-    const scalarField& source,
+    solveScalarField& psi,
+    const solveScalarField& source,
     const direction cmpt
 ) const
 {
@@ -295,6 +296,23 @@ Foam::solverPerformance Foam::petscSolver::solve
 }
 
 
+Foam::solverPerformance Foam::petscSolver::solve
+(
+    scalarField& psi_s,
+    const scalarField& source,
+    const direction cmpt
+) const
+{
+    PrecisionAdaptor<solveScalar, scalar> tpsi(psi_s);
+    return scalarSolve
+    (
+        tpsi.ref(),
+        ConstPrecisionAdaptor<solveScalar, scalar>(source)(),
+        cmpt
+    );
+}
+
+
 void Foam::petscSolver::updateKsp
 (
     KSP& ksp,
diff --git a/src/petsc4Foam/solvers/petscSolver.H b/src/petsc4Foam/solvers/petscSolver.H
index a15f28177e132206501b7ad83f738e079edd15b5..739aac204aa7550a3c2fe33e0f627e5157dc5b8c 100644
--- a/src/petsc4Foam/solvers/petscSolver.H
+++ b/src/petsc4Foam/solvers/petscSolver.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
-    Copyright (C) 2019 Simone Bna
+    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2020 Simone Bna
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -150,12 +150,20 @@ public:
 
     // Member Functions
 
+        //- Solve with given field and rhs (in solveScalar precision).
+        virtual solverPerformance scalarSolve
+        (
+            solveScalarField& psi,
+            const solveScalarField& source,
+            const direction cmpt = 0
+        ) const;
+
         //- Solve the matrix with this solver
         virtual solverPerformance solve
         (
             scalarField& psi,
             const scalarField& source,
-            const direction cmpt=0
+            const direction cmpt = 0
         ) const;
 };
 
diff --git a/src/petsc4Foam/utils/petscWrappedVector.H b/src/petsc4Foam/utils/petscWrappedVector.H
index 97f12124f0b66a009690bf6faea07be868aff021..88fd2dc0cb618e091ad8e00e9c93424d853691e7 100644
--- a/src/petsc4Foam/utils/petscWrappedVector.H
+++ b/src/petsc4Foam/utils/petscWrappedVector.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -78,7 +78,7 @@ class PetscWrappedVector
             VecPlaceArray
             (
                 vec_,
-                reinterpret_cast<scalar*>(const_cast<T*>(list.cdata()))
+                reinterpret_cast<solveScalar*>(const_cast<T*>(list.cdata()))
             );
             reset_ = true;
         }
@@ -97,7 +97,7 @@ class PetscWrappedVector
                 pTraits<T>::nComponents,
                 list.size(),
                 PETSC_DECIDE,
-                reinterpret_cast<scalar*>(const_cast<T*>(list.cdata())),
+                reinterpret_cast<solveScalar*>(const_cast<T*>(list.cdata())),
                 &vec_
             );
             reset_ = false;
@@ -130,7 +130,7 @@ public:
         //- Wrap an OpenFOAM list of values as a PETSc vector
         PetscWrappedVector
         (
-            const UList<scalar>& list,
+            const UList<solveScalar>& list,
             Mat mat
         )
         :
@@ -142,7 +142,7 @@ public:
         //- Wrap an OpenFOAM list of values as a PETSc vector
         explicit PetscWrappedVector
         (
-            const UList<scalar>& list,
+            const UList<solveScalar>& list,
             MPI_Comm comm = PETSC_COMM_WORLD
         )
         :