diff --git a/CMakeLists.txt b/CMakeLists.txt
index d8dcbd8e3..720a80fdf 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -87,6 +87,9 @@ elseif (OGS_CONFIG STREQUAL GEMS)
 	set(OGS_CONFIG_DEFAULT_CHEMSOLVER "GEMS" CACHE INTERNAL "" FORCE)
 elseif (OGS_CONFIG STREQUAL IPQC)
 	set(OGS_CONFIG_DEFAULT_CHEMSOLVER "IPQC" CACHE INTERNAL "" FORCE)
+	if(PARALLEL_USE_MPI)
+		set(OGS_CONFIG_DEFAULT_LSOLVER "SP" CACHE INTERNAL "" FORCE)
+	endif()
 elseif (OGS_CONFIG STREQUAL PQC)
 	set(OGS_CONFIG_DEFAULT_LSOLVER "SP" CACHE INTERNAL "" FORCE)
 	set(OGS_CONFIG_DEFAULT_CHEMSOLVER "PQC" CACHE INTERNAL "" FORCE)
diff --git a/FEM/fem_ele_std.cpp b/FEM/fem_ele_std.cpp
index a34211dd1..14ae5dfaa 100644
--- a/FEM/fem_ele_std.cpp
+++ b/FEM/fem_ele_std.cpp
@@ -1578,9 +1578,9 @@ double CFiniteElementStd::CalCoefMass()
             val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta);
             double drho_dp_rho = 0.0;
             // get drho/dp/rho from material model or direct input
+            const double rho_val = FluidProp->Density();
             if (FluidProp->compressibility_model_pressure > 0)
             {
-                const double rho_val = FluidProp->Density();
                 double arg[2];
                 arg[0] = interpolate(NodalVal1);   //   p
                 arg[1] = interpolate(NodalValC1);  //   T
@@ -1588,7 +1588,10 @@ double CFiniteElementStd::CalCoefMass()
             }
             else
             {
-                drho_dp_rho = FluidProp->drho_dp;
+                drho_dp_rho =
+                    rho_val > DBL_EPSILON
+                        ? FluidProp->drho_dp * FluidProp->rho_0 / rho_val
+                        : 0.0;
             }
 
             const double poro_val =
@@ -1603,7 +1606,7 @@ double CFiniteElementStd::CalCoefMass()
                 h_fem = this;
                 storage_effstress =
                     MediaProp->StorageFunctionEffStress(Index, nnodes, h_fem);
-                val *= storage_effstress;
+                val += storage_effstress;
             }
 
             // val /= time_unit_factor;
diff --git a/FEM/rf_mmp_new.cpp b/FEM/rf_mmp_new.cpp
index 6a85729e7..814f05242 100644
--- a/FEM/rf_mmp_new.cpp
+++ b/FEM/rf_mmp_new.cpp
@@ -8001,7 +8001,7 @@ double CMediumProperties::StorageFunction(long index, double* gp, double theta)
     return storage;
 }
 
-// AS:08.2012 storgae function eff stress
+// AS:08.2012 storage function of effective stress
 double CMediumProperties::StorageFunctionEffStress(long index, int nnodes,
                                                    CFiniteElementStd* h_fem)
 {