From f2b68d8c34fa02a3e4be21b7cd24e7c4f2f17d7c Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Thu, 14 Apr 2016 15:21:54 -0400
Subject: [PATCH] [Examples] General cleanup of C++ flamespeed example

---
 samples/cxx/flamespeed/flamespeed.cpp         | 113 +++++++-----------
 .../cxx/flamespeed/flamespeed_blessed_0.csv   |   3 -
 2 files changed, 45 insertions(+), 71 deletions(-)

diff --git a/samples/cxx/flamespeed/flamespeed.cpp b/samples/cxx/flamespeed/flamespeed.cpp
index b23fc8c25..f612d4b86 100644
--- a/samples/cxx/flamespeed/flamespeed.cpp
+++ b/samples/cxx/flamespeed/flamespeed.cpp
@@ -8,10 +8,10 @@
 #include "cantera/oneD/StFlow.h"
 #include "cantera/IdealGasMix.h"
 #include "cantera/transport.h"
+#include <fstream>
 
 using namespace Cantera;
-using std::cout;
-using std::endl;
+using fmt::print;
 
 int flamespeed(double phi)
 {
@@ -20,45 +20,31 @@ int flamespeed(double phi)
 
         doublereal temp = 300.0; // K
         doublereal pressure = 1.0*OneAtm; //atm
-        doublereal uin=0.3; //m/sec
+        doublereal uin = 0.3; //m/sec
 
-        gas.setState_TPX(temp, pressure, "CH4:1.0, O2:2.0, N2:7.52");
         size_t nsp = gas.nSpecies();
+        vector_fp x(nsp, 0.0);
 
-        vector_fp x(nsp);
+        doublereal C_atoms = 1.0;
+        doublereal H_atoms = 4.0;
+        doublereal ax = C_atoms + H_atoms / 4.0;
+        doublereal fa_stoic = 1.0 / (4.76 * ax);
+        x[gas.speciesIndex("CH4")] = 1.0;
+        x[gas.speciesIndex("O2")] = 0.21 / phi / fa_stoic;
+        x[gas.speciesIndex("N2")] = 0.79 / phi/ fa_stoic;
 
-        doublereal C_atoms=1.0;
-        doublereal H_atoms=4.0;
-        doublereal ax=C_atoms+H_atoms/4.0;
-        doublereal fa_stoic=1.0/(4.76*ax);
-        for (size_t k=0; k<nsp; k++) {
-            if (k==gas.speciesIndex("CH4")) {
-                x[k]=1.0;
-            } else if (k==gas.speciesIndex("O2")) {
-                x[k]=0.21/phi/fa_stoic;
-            } else if (k==gas.speciesIndex("N2")) {
-                x[k]=0.79/phi/fa_stoic;
-            } else {
-                x[k]=0.0;
-            }
-        }
-
-        gas.setState_TPX(temp,pressure,x.data());
-        doublereal rho_in=gas.density();
+        gas.setState_TPX(temp, pressure, x.data());
+        doublereal rho_in = gas.density();
 
         vector_fp yin(nsp);
         gas.getMassFractions(&yin[0]);
 
-        try {
-            gas.equilibrate("HP");
-        } catch (CanteraError& err) {
-            std::cout << err.what() << std::endl;
-        }
+        gas.equilibrate("HP");
         vector_fp yout(nsp);
         gas.getMassFractions(&yout[0]);
         doublereal rho_out = gas.density();
-        doublereal Tad=gas.temperature();
-        cout << phi<<' '<<Tad<<endl;
+        doublereal Tad = gas.temperature();
+        print("phi = {}, Tad = {}\n", phi, Tad);
 
         //=============  build each domain ========================
 
@@ -69,11 +55,11 @@ int flamespeed(double phi)
 
         // create an initial grid
         int nz = 6;
-        doublereal lz=0.1;
+        doublereal lz = 0.1;
         vector_fp z(nz);
-        doublereal dz=lz/((doublereal)(nz-1));
-        for (int iz=0; iz<nz; iz++) {
-            z[iz]=((doublereal)iz)*dz;
+        doublereal dz = lz/((doublereal)(nz-1));
+        for (int iz = 0; iz < nz; iz++) {
+            z[iz] = ((doublereal)iz)*dz;
         }
 
         flow.setupGrid(nz, &z[0]);
@@ -112,7 +98,7 @@ int flamespeed(double phi)
         vector_fp locs{0.0, 0.3, 0.7, 1.0};
         vector_fp value;
 
-        double uout=inlet.mdot()/rho_out;
+        double uout = inlet.mdot()/rho_out;
         value = {uin, uin, uout, uout};
         flame.setInitialGuess("u",locs,value);
         value = {temp, temp, Tad, Tad};
@@ -129,10 +115,10 @@ int flamespeed(double phi)
 
         flame.showSolution();
 
-        int flowdomain=1;
-        double ratio=10.0;
-        double slope=0.08;
-        double curve=0.1;
+        int flowdomain = 1;
+        double ratio = 10.0;
+        double slope = 0.08;
+        double curve = 0.1;
 
         flame.setRefineCriteria(flowdomain,ratio,slope,curve);
 
@@ -149,58 +135,49 @@ int flamespeed(double phi)
 
         flame.solve(loglevel,refine_grid);
         double flameSpeed_mix = flame.value(flowdomain,flow.componentIndex("u"),0);
-        cout << "Flame speed with mixture-averaged transport: " <<
-             flame.value(flowdomain,flow.componentIndex("u"),0) << " m/s" << endl;
+        print("Flame speed with mixture-averaged transport: {} m/s\n",
+              flameSpeed_mix);
 
         // now switch to multicomponent transport
         flow.setTransport(*trmulti);
         flame.solve(loglevel, refine_grid);
         double flameSpeed_multi = flame.value(flowdomain,flow.componentIndex("u"),0);
-        cout << "Flame speed with multicomponent transport: " <<
-             flame.value(flowdomain,flow.componentIndex("u"),0) << " m/s" << endl;
+        print("Flame speed with multicomponent transport: {} m/s\n",
+              flameSpeed_multi);
 
         // now enable Soret diffusion
         flow.enableSoret(true);
         flame.solve(loglevel, refine_grid);
         double flameSpeed_full = flame.value(flowdomain,flow.componentIndex("u"),0);
-        cout << "Flame speed with multicomponent transport + Soret: " <<
-             flame.value(flowdomain,flow.componentIndex("u"),0) << " m/s" << endl;
+        print("Flame speed with multicomponent transport + Soret: {} m/s\n",
+              flameSpeed_full);
 
-        int np=flow.nPoints();
         vector_fp zvec,Tvec,COvec,CO2vec,Uvec;
 
-        printf("\n%9s\t%8s\t%5s\t%7s\n","z (m)", "T (K)", "U (m/s)", "Y(CO)");
-        for (int n=0; n<np; n++) {
+        print("\n{:9s}\t{:8s}\t{:5s}\t{:7s}\n",
+              "z (m)", "T (K)", "U (m/s)", "Y(CO)");
+        for (size_t n = 0; n < flow.nPoints(); n++) {
             Tvec.push_back(flame.value(flowdomain,flow.componentIndex("T"),n));
             COvec.push_back(flame.value(flowdomain,flow.componentIndex("CO"),n));
             CO2vec.push_back(flame.value(flowdomain,flow.componentIndex("CO2"),n));
             Uvec.push_back(flame.value(flowdomain,flow.componentIndex("u"),n));
             zvec.push_back(flow.grid(n));
-            printf("%9.6f\t%8.3f\t%5.3f\t%7.5f\n",flow.grid(n),Tvec[n],Uvec[n],COvec[n]);
+            print("{:9.6f}\t{:8.3f}\t{:5.3f}\t{:7.5f}\n",
+                  flow.grid(n), Tvec[n], Uvec[n], COvec[n]);
         }
 
-        cout << endl<<"Adiabatic flame temperature from equilibrium is: "<<Tad<<endl;
-        cout << "Flame speed for phi="<<phi<<" is "<<Uvec[0]<<" m/s."<<endl;
+        print("\nAdiabatic flame temperature from equilibrium is: {}\n", Tad);
+        print("Flame speed for phi={} is {} m/s.\n", phi, Uvec[0]);
 
-        std::string reportFile = "flamespeed.csv";
-        FILE* FP = fopen(reportFile.c_str(), "w");
-        if (!FP) {
-            printf("Failure to open file\n");
-            exit(-1);
+        std::ofstream outfile("flamespeed.csv", std::ios::trunc);
+        outfile << "  Grid,   Temperature,   Uvec,   CO,    CO2\n";
+        for (size_t n = 0; n < flow.nPoints(); n++) {
+            print(outfile, " {:11.3e}, {:11.3e}, {:11.3e}, {:11.3e}, {:11.3e}\n",
+                  flow.grid(n), Tvec[n], Uvec[n], COvec[n], CO2vec[n]);
         }
-
-        fprintf(FP," Flame speed (mixture-averaged      ) = %11.3e m/s\n", flameSpeed_mix);
-        fprintf(FP," Flame speed (multicomponent        ) = %11.3e m/s\n", flameSpeed_multi);
-        fprintf(FP," Flame speed (multicomponent + Soret) = %11.3e m/s\n", flameSpeed_full);
-        fprintf(FP,"  Grid,   Temperature,   Uvec,   CO,    CO2\n");
-        for (int n = 0; n < np; n++) {
-            fprintf(FP," %11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
-                    flow.grid(n), Tvec[n], Uvec[n], COvec[n], CO2vec[n]);
-        }
-        fclose(FP);
     } catch (CanteraError& err) {
         std::cerr << err.what() << std::endl;
-        std::cerr << "program terminating." << endl;
+        std::cerr << "program terminating." << std::endl;
         return -1;
     }
     return 0;
@@ -209,7 +186,7 @@ int flamespeed(double phi)
 int main()
 {
     double phi;
-    cout << "Enter phi: ";
+    print("Enter phi: ");
     std::cin >> phi;
     return flamespeed(phi);
 }
diff --git a/samples/cxx/flamespeed/flamespeed_blessed_0.csv b/samples/cxx/flamespeed/flamespeed_blessed_0.csv
index 370816ad8..4f6b9bcde 100644
--- a/samples/cxx/flamespeed/flamespeed_blessed_0.csv
+++ b/samples/cxx/flamespeed/flamespeed_blessed_0.csv
@@ -1,6 +1,3 @@
- Flame speed (mixture-averaged      ) =   3.461e-01 m/s
- Flame speed (multicomponent        ) =   3.516e-01 m/s
- Flame speed (multicomponent + Soret) =   3.516e-01 m/s
   Grid,   Temperature,   Uvec,   CO,    CO2
    0.000e+00,   3.000e+02,   3.516e-01,   4.128e-11,   2.311e-12
    1.202e-03,   3.000e+02,   3.516e-01,   8.802e-10,   6.402e-11