diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
index 16ea5af785cd1e0ea8613d841e6f227d8f9dec66..3cac582cf307b26322de085ef3383b4050dcd5ed 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H
@@ -253,6 +253,7 @@ public:
             //- Read and reset the solver parameters from the given stream
             virtual void read(const dictionary&);
 
+            //- Solve with given field and rhs
             virtual solverPerformance solve
             (
                 scalarField& psi,
@@ -260,6 +261,15 @@ public:
                 const direction cmpt=0
             ) const = 0;
 
+            //- Solve with given field and rhs (in solveScalar precision).
+            //  Default is to call solve routine
+            virtual solverPerformance scalarSolve
+            (
+                solveScalarField& psi,
+                const solveScalarField& source,
+                const direction cmpt=0
+            ) const;
+
             //- Return the matrix norm used to normalise the residual for the
             //- stopping criterion
             solveScalarField::cmptType normFactor
diff --git a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
index 42e24d04b5b1b36ad7d4e3fa1cbb42101f6070b8..50b6f46bcc8638d095a78fb52768f6c0140dad55 100644
--- a/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
+++ b/src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C
@@ -27,6 +27,7 @@ License
 
 #include "lduMatrix.H"
 #include "diagonalSolver.H"
+#include "PrecisionAdaptor.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -170,6 +171,23 @@ void Foam::lduMatrix::solver::read(const dictionary& solverControls)
 }
 
 
+Foam::solverPerformance Foam::lduMatrix::solver::scalarSolve
+(
+    solveScalarField& psi,
+    const solveScalarField& source,
+    const direction cmpt
+) const
+{
+    PrecisionAdaptor<scalar, solveScalar> tpsi_s(psi);
+    return solve
+    (
+        tpsi_s.ref(),
+        ConstPrecisionAdaptor<scalar, solveScalar>(source)(),
+        cmpt
+    );
+}
+
+
 Foam::solveScalarField::cmptType Foam::lduMatrix::solver::normFactor
 (
     const solveScalarField& psi,
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
index 644de4eef9a94fad2f978991fb83af796d83e03b..d5cbc20296eb39393e3de143a44bc751acdb53cf 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C
@@ -695,7 +695,7 @@ void Foam::GAMGSolver::solveCoarsestLevel
         coarsestCorrField = 0;
         const solverPerformance coarseSolverPerf
         (
-            coarsestSolverPtr_->solve
+            coarsestSolverPtr_->scalarSolve
             (
                 coarsestCorrField,
                 coarsestSource