diff --git a/test/material/pengrobinson/test_pengrobinson.cc b/test/material/pengrobinson/test_pengrobinson.cc
index 508b566b8e910a6f74c6b04928f8e603aba91979..3996032d46ff998b98edb79179253504692a16c7 100644
--- a/test/material/pengrobinson/test_pengrobinson.cc
+++ b/test/material/pengrobinson/test_pengrobinson.cc
@@ -23,6 +23,7 @@
  *        Peng-Robinson EOS) and the NCP flash solver.
  */
 #include <config.h>
+#include <iomanip>
 
 #include <dumux/material/constraintsolvers/ncpflash.hh>
 #include <dumux/material/fluidstates/compositional.hh>
@@ -217,7 +218,7 @@ int main(int argc, char** argv)
     Scalar rho_gRef = surfaceFluidState.density(gPhaseIdx);
     Scalar rho_oRef = surfaceFluidState.density(oPhaseIdx);
 
-    std::cout << "alpha[-] p[Pa] S_g[-] rho_o[kg/m^3] rho_g[kg/m^3] <M_o>[kg/mol] <M_g>[kg/mol] R_s[m^3/m^3] B_g[-] B_o[-]\n";
+    std::cout << "alpha[-]     p[Pa]        S_g[-]       rho_o[kg/m^3] rho_g[kg/m^3] rhom_o[mol/m^3] rhom_g[mol/m^3] <M_o>[kg/mol] <M_g>[kg/mol] R_s[m^3/m^3] B_g[-] B_o[-]\n";
     int n = 50;
     for (int i = 0; i < n; ++i) {
         Scalar minAlpha = 0.98;
@@ -237,11 +238,13 @@ int main(int argc, char** argv)
                                                                      surfaceAlpha,
                                                                      flashFluidState,
                                                                      /*guessInitial=*/false);
-        std::cout << alpha << " "
+        std::cout << std::scientific << alpha << " "
                   << flashFluidState.pressure(oPhaseIdx) << " "
                   << flashFluidState.saturation(gPhaseIdx) << " "
                   << flashFluidState.density(oPhaseIdx) << " "
                   << flashFluidState.density(gPhaseIdx) << " "
+                  << flashFluidState.molarDensity(oPhaseIdx) << " "
+                  << flashFluidState.molarDensity(gPhaseIdx) << " "
                   << flashFluidState.averageMolarMass(oPhaseIdx) << " "
                   << flashFluidState.averageMolarMass(gPhaseIdx) << " "
                   << surfaceFluidState.saturation(gPhaseIdx)*alphaSurface << " "