diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
index 65143f5d27..84ba7e84cd 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C
@@ -232,12 +232,22 @@ void reactingOneDim::solveContinuity()
         Info<< "reactingOneDim::solveContinuity()" << endl;
     }
 
-    solve
-    (
-        fvm::ddt(rho_)
-     ==
-      - solidChemistry_->RRg()
-    );
+    if (moveMesh_)
+    {
+        const scalarField mass0 = rho_*regionMesh().V();
+
+        fvScalarMatrix rhoEqn
+        (
+            fvm::ddt(rho_)
+         ==
+          - solidChemistry_->RRg()
+        );
+
+        rhoEqn.solve();
+
+        updateMesh(mass0);
+
+    }
 }
 
 
@@ -261,14 +271,15 @@ void reactingOneDim::solveSpeciesMass()
             solidChemistry_->RRs(i)
         );
 
-        if (moveMesh_)
+        if (regionMesh().moving())
         {
-            surfaceScalarField phiRhoMesh
+            surfaceScalarField phiYiRhoMesh
             (
                 fvc::interpolate(Yi*rho_)*regionMesh().phi()
             );
 
-            YiEqn -= fvc::div(phiRhoMesh);
+            YiEqn += fvc::div(phiYiRhoMesh);
+
         }
 
         YiEqn.solve(regionMesh().solver("Yi"));
@@ -303,14 +314,14 @@ void reactingOneDim::solveEnergy()
       + fvc::div(phiGas)
     );
 
-    if (moveMesh_)
+    if (regionMesh().moving())
     {
-        surfaceScalarField phiMesh
+        surfaceScalarField phihMesh
         (
             fvc::interpolate(rho_*h_)*regionMesh().phi()
         );
 
-        hEqn -= fvc::div(phiMesh);
+        hEqn += fvc::div(phihMesh);
     }
 
     hEqn.relax();
@@ -349,7 +360,18 @@ reactingOneDim::reactingOneDim(const word& modelType, const fvMesh& mesh)
     pyrolysisModel(modelType, mesh),
     solidChemistry_(solidChemistryModel::New(regionMesh())),
     solidThermo_(solidChemistry_->solid()),
-    rho_(solidThermo_.rho()),
+    rho_
+    (
+        IOobject
+        (
+            "rho",
+            regionMesh().time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        solidThermo_.rho()
+    ),
     Ys_(solidThermo_.composition().Y()),
     h_(solidThermo_.he()),
     primaryRadFluxName_(coeffs().lookupOrDefault<word>("radFluxName", "Qr")),
@@ -449,7 +471,18 @@ reactingOneDim::reactingOneDim
     pyrolysisModel(modelType, mesh, dict),
     solidChemistry_(solidChemistryModel::New(regionMesh())),
     solidThermo_(solidChemistry_->solid()),
-    rho_(solidThermo_.rho()),
+    rho_
+    (
+        IOobject
+        (
+            "rho",
+            regionMesh().time().timeName(),
+            regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        solidThermo_.rho()
+    ),
     Ys_(solidThermo_.composition().Y()),
     h_(solidThermo_.he()),
     primaryRadFluxName_(dict.lookupOrDefault<word>("radFluxName", "Qr")),
@@ -681,17 +714,15 @@ void reactingOneDim::evolveRegion()
 {
     Info<< "\nEvolving pyrolysis in region: " << regionMesh().name() << endl;
 
-    const scalarField mass0 = rho_*regionMesh().V();
-
     solidChemistry_->solve
     (
         time().value() - time().deltaTValue(),
         time().deltaTValue()
     );
 
-    solveContinuity();
+    calculateMassTransfer();
 
-    updateMesh(mass0);
+    solveContinuity();
 
     chemistrySh_ = solidChemistry_->Sh()();
 
@@ -704,9 +735,9 @@ void reactingOneDim::evolveRegion()
         solveEnergy();
     }
 
-    calculateMassTransfer();
-
     solidThermo_.correct();
+
+    rho_ = solidThermo_.rho();
 }
 
 
diff --git a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
index aa3d7baa9b..619c9c0e9a 100644
--- a/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
+++ b/src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H
@@ -84,14 +84,8 @@ protected:
 
         // Reference to solid thermo properties
 
-//             //- Absorption coefficient [1/m]
-//             const volScalarField& kappaRad_;
-//
-//             //- Thermal conductivity [W/m/K]
-//             const volScalarField& kappa_;
-
             //- Density [kg/m3]
-            volScalarField& rho_;
+            volScalarField rho_;
 
             //- List of solid components
             PtrList<volScalarField>& Ys_;
@@ -221,8 +215,8 @@ public:
 
             //- Fields
 
-                //- Return density [kg/m3]
-                virtual const volScalarField& rho() const;
+                //- Return const density [Kg/m3]
+                const volScalarField& rho() const;
 
                 //- Return const temperature [K]
                 virtual const volScalarField& T() const;
diff --git a/src/regionModels/regionModel/regionModel1D/regionModel1D.C b/src/regionModels/regionModel/regionModel1D/regionModel1D.C
index 372e19a3f3..60243b55cd 100644
--- a/src/regionModels/regionModel/regionModel1D/regionModel1D.C
+++ b/src/regionModels/regionModel/regionModel1D/regionModel1D.C
@@ -117,6 +117,8 @@ void Foam::regionModels::regionModel1D::initialise()
     }
 
     boundaryFaceOppositeFace_.setSize(localPyrolysisFaceI);
+    boundaryFaceFaces_.setSize(localPyrolysisFaceI);
+    boundaryFaceCells_.setSize(localPyrolysisFaceI);
 
     surfaceScalarField& nMagSf = nMagSfPtr_();
 
@@ -192,6 +194,7 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 
     const polyBoundaryMesh& bm = regionMesh().boundaryMesh();
 
+    label totalFaceId = 0;
     forAll(intCoupledPatchIDs_, localPatchI)
     {
         label patchI = intCoupledPatchIDs_[localPatchI];
@@ -200,8 +203,9 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 
         forAll(pp, patchFaceI)
         {
-            const labelList& faces = boundaryFaceFaces_[patchFaceI];
-            const labelList& cells = boundaryFaceCells_[patchFaceI];
+            const labelList& faces = boundaryFaceFaces_[totalFaceId];
+            const labelList& cells = boundaryFaceCells_[totalFaceId];
+
             const vector n = pp.faceNormals()[patchFaceI];
             const vector sf = pp.faceAreas()[patchFaceI];
 
@@ -231,7 +235,7 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
 
                     if
                     (
-                        ((nbrCf - (oldPoints[pointI] + newDelta)) & n)
+                        mag((nbrCf - (oldPoints[pointI] + newDelta)) & n)
                       > minDelta
                     )
                     {
@@ -242,19 +246,17 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
                 }
                 nbrCf = oldCf[i + 1] + localDelta;
             }
-
             // Modify boundary
-            const label bFaceI = boundaryFaceOppositeFace_[patchFaceI];
+            const label bFaceI = boundaryFaceOppositeFace_[totalFaceId];
             const face f = regionMesh().faces()[bFaceI];
             const label cellI = cells[cells.size() - 1];
             newDelta += (deltaV[cellI]/mag(sf))*n;
             forAll(f, pti)
             {
                 const label pointI = f[pti];
-
                 if
                 (
-                    ((nbrCf - (oldPoints[pointI] + newDelta)) & n)
+                    mag((nbrCf - (oldPoints[pointI] + newDelta)) & n)
                   > minDelta
                 )
                 {
@@ -262,9 +264,9 @@ Foam::tmp<Foam::labelField> Foam::regionModels::regionModel1D::moveMesh
                     cellMoveMap[cellI] = 1;
                 }
             }
+            totalFaceId ++;
         }
     }
-
     // Move points
     regionMesh().movePoints(newPoints);