//#r "D:\BoSSS2\experimental\public\src\L4-application\BoSSSpad\bin\Release\net5.0\bossspad.dll"
#r "BoSSSpad.dll"
using System;
using System.Collections.Generic;
using System.Linq;
using System.IO;
using System.Data;
using System.Globalization;
using System.Threading;
using ilPSP;
using ilPSP.Utils;
using BoSSS.Platform;
using BoSSS.Foundation;
using BoSSS.Foundation.Grid;
using BoSSS.Foundation.Grid.Classic;
using BoSSS.Foundation.IO;
using BoSSS.Solution;
using BoSSS.Solution.Control;
using BoSSS.Solution.GridImport;
using BoSSS.Solution.Statistic;
using BoSSS.Solution.Utils;
using BoSSS.Solution.Gnuplot;
using BoSSS.Application.BoSSSpad;
using BoSSS.Application.XNSE_Solver;
using static BoSSS.Application.BoSSSpad.BoSSSshell;
using BoSSS.Foundation.Grid.RefElements;
using BoSSS.Platform.LinAlg;
using BoSSS.Solution.NSECommon;
using BoSSS.Application.XNSEC;
Init();

The below script needs to be able to find the current output cell; this is an easy method to get it.

Open Database

BoSSSshell.WorkflowMgm.Init("CounterFlowFlame_MF_FullComparison");

Project name is set to 'CounterFlowFlame_MF_FullComparison'.

Default Execution queue is chosen for the database.

Opening existing database '\\fdygitrunner\ValidationTests\databases\CounterFlowFlame_MF_FullComparison'.

static var myBatch = GetDefaultQueue();
static var myDb = BoSSSshell.WorkflowMgm.DefaultDatabase;

BoSSSshell.WorkflowMgm.SetNameBasedSessionJobControlCorrelation()

Case configuration

int nCells = 10;
int[] multiplierS = new int[]{2,5,11}; // multiplying factor of inlet velocities
// int[] multiplierS = new int[]{2}; // multiplying factor of inlet velocities
double[] Smoothings = new double[] {80};
int dgMF = 3; // Mass fraction DG degree
int[] DGDegrees = new int[] {4};

int[] nCellsArray = new int[] {10};

double[] xCellsMultiplier = new double[] {3};

double InitialMassFuelIn = 0.02400; //kg/m2s 
double InitialMassAirIn = 0.02400 *3; //kg/m2s 
bool parabolicVelocityProfile = false;
bool chemicalReactionActive = true; 
bool useFullGeometry = false;
bool wallBounded = true;

Create grid

public static class GridFactory {

    public static Grid2D GenerateGrid(int nCells, bool wallBounded, double xDensity) {
        double L = 0.02;
        double LRef = L;
        double R_dim = L / 2;
        double separation = 1.0; // nondimensional
        double xleft = 0;
        double xright = separation;
        double radius_inlet = R_dim / LRef;
        double ybot = -separation * 3; 
        double ytop = separation * 3; 

        // Define Grid
        double[] _xNodes;
        double[] _yNodes;

       // X-NODES
        double xNodesDensity = xDensity; 
        _xNodes = GenericBlas.Linspace(xleft, xright, (int)(nCells * xNodesDensity + 1));
        double sf3 = 0.80;
        var myCenterNodes = GenericBlas.SinLinSpacing(0, radius_inlet*2,0.7, nCells*2 + 1).ToList(); // center
        var centerNodes = myCenterNodes.GetSubVector(0, myCenterNodes.Count / 2 + 1); // Take only "bottom side" of node array


        List<double> yNodesTop = (GenericBlas.SinLinSpacing(radius_inlet, (ytop - radius_inlet) * 2 + radius_inlet, sf3, nCells * 4 + 1).ToList()); // Nodes corresponding to the oxidizer inlet, right part
        var myYnodesTop = yNodesTop.GetSubVector(0, yNodesTop.Count / 2 + 1); // Take only "bottom side" of node array

        List<double> yUpperPart = new List<double>();
        yUpperPart.AddRange(centerNodes.SkipLast(1));
        yUpperPart.AddRange(myYnodesTop.SkipLast(0));

        var yBottomPart = (yUpperPart.ToArray()).CloneAs();
        yBottomPart.ScaleV(-1.0);
        Array.Reverse(yBottomPart);
        var yNodes = new List<double>();
        yNodes.AddRange(yBottomPart.SkipLast(1));
        yNodes.AddRange(yUpperPart.SkipLast(0));
        
        var grd = Grid2D.Cartesian2DGrid(_xNodes, yNodes.ToArray());
        // Define Edge tags
        grd.EdgeTagNames.Add(1, "Velocity_Inlet_CH4");
        grd.EdgeTagNames.Add(2, "Velocity_Inlet_O2");
        grd.EdgeTagNames.Add(3, "Pressure_Outlet");

        if(wallBounded){
            grd.EdgeTagNames.Add(4, "Wall");
        }
        grd.DefineEdgeTags(delegate (double[] X) {
            double x = X[0];
            double y = X[1];
            //Edge tags
            //1: Velocity inlet O_2
            //2: Velocity inlet CH_4
            //3: Pressure outlet
            if(Math.Abs(x - xleft) < 1e-8) { // Left boundary
                if(Math.Abs(y) - radius_inlet < 1e-8) { // Fuel Inlet
                    return 1;
                } else {
                    if(wallBounded){
                    return 4;
                } else{
                    return 3;
                }
                }
            }
            if(Math.Abs(x - xright) < 1e-8) { // right boundary
                if(Math.Abs(y) - radius_inlet < 1e-8) { // oxy Inlet
                    return 2;//2
                } else {
                    if(wallBounded){
                        return 4;
                    } else{
                        return 3;
                    }
                }
            } else {
                return 3; //3 Pressure outlet
            }
        });
        myDb.SaveGrid(ref grd);

        return grd;
    }

    public static Grid2D GenerateGrid_FullGeometry(int nCells) {
        double L = 0.02;
        double LRef = L;
        double R_dim = L / 2;
        double separation = 1.0; // nondimensional
        double xleft = 0;
        double xright = separation;
        double radius_inlet = R_dim / LRef;
        double ybot = -separation * 3; //
        double ytop = separation * 3;//

        // Define Grid
        double[] _xNodes;
        double[] _yNodes;

         // X-NODES
        double[] XNODES1 = GenericBlas.Linspace(0, separation, (int)(nCells / 2 + 1));
        double[] XNODES2 = GenericBlas.Linspace(separation * 1, separation * 2, (int)(nCells + 1));
        double[] XNODES3 = GenericBlas.Linspace(separation * 2, separation * 3, (int)(nCells / 2 + 1));
        List<double> xNodesList = new List<double>();
        xNodesList.AddRange(XNODES1.Take(XNODES1.Length - 1));
        xNodesList.AddRange(XNODES2.Take(XNODES2.Length - 1));
        xNodesList.AddRange(XNODES3);
        _xNodes = xNodesList.ToArray();

        // Y-NODES
        double sf3 = 0.0;
        var _yNodes2 = GenericBlas.Linspace(-radius_inlet, radius_inlet, nCells + 1).ToList(); // center
        List<double> yNodesTop = (GenericBlas.SinLinSpacing(radius_inlet, (ytop - radius_inlet) * 2 + radius_inlet, sf3, nCells * 2 + 1).ToList()); // Nodes corresponding to the oxidizer inlet, right part
        var myYnodesTop = yNodesTop.GetSubVector(0, yNodesTop.Count / 2 + 1); // Take only "bottom side" of node array
        var myYNodesBot = myYnodesTop.CloneAs();
        myYNodesBot.ScaleV(-1.0);
        Array.Reverse(myYNodesBot);
        List<double> list2 = new List<double>();
        list2.AddRange(myYNodesBot.Take(myYNodesBot.Length - 1).ToList());
        list2.AddRange(_yNodes2.Take(_yNodes2.Count() - 1).ToList());
        list2.AddRange(myYnodesTop.ToList());
        _yNodes = list2.ToArray();

        double[] CutOut1Point1 = new double[2] { separation, separation / 2 };
        double[] CutOut1Point2 = new double[2] { 0, -separation / 2 };

        double[] CutOut2Point1 = new double[2] { 2 * separation, -separation / 2 };
        double[] CutOut2Point2 = new double[2] { 3 * separation, separation / 2 };

        var CutOut1 = new BoSSS.Platform.Utils.Geom.BoundingBox(2);
        CutOut1.AddPoint(CutOut1Point1);
        CutOut1.AddPoint(CutOut1Point2);

        var CutOut2 = new BoSSS.Platform.Utils.Geom.BoundingBox(2);
        CutOut2.AddPoint(CutOut2Point1);
        CutOut2.AddPoint(CutOut2Point2);
        var CutOuts = new BoSSS.Platform.Utils.Geom.BoundingBox[] { CutOut1, CutOut2 };
        var grd = Grid2D.Cartesian2DGrid(_xNodes, _yNodes, CutOuts: CutOuts);

        grd.EdgeTagNames.Add(1, "Velocity_Inlet_CH4");
        grd.EdgeTagNames.Add(2, "Velocity_Inlet_O2");
        grd.EdgeTagNames.Add(3, "Pressure_Outlet");
        grd.EdgeTagNames.Add(4, "Wall");

        grd.DefineEdgeTags(delegate (double[] X) {
            double x = X[0];
            double y = X[1];
            //Edge tags
            //1: Velocity inlet O_2
            //2: Velocity inlet CH_4
            //3: Pressure outlet
            if(Math.Abs(x - separation) < 1e-8 && Math.Abs(y) - radius_inlet < 1e-8) {
                return 1;
            }
            if(Math.Abs(x - 2 * separation) < 1e-8 && Math.Abs(y) - radius_inlet < 1e-8) {
                return 2;
            }
            if(Math.Abs(y - radius_inlet) < 1e-8 && x > 0 && x < separation) {
                return 4;
            }
            if(Math.Abs(y - radius_inlet) < 1e-8 && x > 2 * separation && x < 3 * separation) {
                return 4;
            }
            if(Math.Abs(y + radius_inlet) < 1e-8 && x > 0 && x < separation) {
                return 4;
            }
            if(Math.Abs(y + radius_inlet) < 1e-8 && x > 2 * separation && x < 3 * separation) {
                return 4;
            }
            return 3;
        }
        );
        myDb.SaveGrid(ref grd);

        return grd;
    }
}

(16,18): warning CS0168: The variable '_yNodes' is declared but never used



(91,16): warning CS0219: The variable 'xleft' is assigned but its value is never used

public static class BoundaryValueFactory { 

    public static string GetPrefixCode(double ConstVal, double inletRadius, double uInFuel, double uInAir, double sigma) {
        using(var stw = new System.IO.StringWriter()) {

            stw.WriteLine("static class BoundaryValues {");

            stw.WriteLine("  static public double ConstantValue(double[] X) {");
            stw.WriteLine("    return "+ ConstVal +";");
            stw.WriteLine("  }");

            stw.WriteLine("  static public double ParabolaVelocityFuel(double[] X) {");
            stw.WriteLine("    return  (1.0 - Math.Pow(X[1] / "+inletRadius+", 2)) * "+uInFuel+" ;");
            stw.WriteLine("  }");

            stw.WriteLine("  static public double ParabolaVelocityAir(double[] X) {");
            stw.WriteLine("    return  -(1.0 - Math.Pow(X[1] / "+inletRadius+", 2)) * "+uInAir+";");
            stw.WriteLine("  }");

            stw.WriteLine("  static public double RegularizedPlugFlowFuel(double[] X) {");    
            stw.WriteLine(" double res = 0;");  
            stw.WriteLine("         if(X[1] > 0) { ");
            stw.WriteLine("     double H = 0.5 * (1.0 + Math.Tanh("+sigma+" * (X[1] - "+inletRadius+"))); ");
            stw.WriteLine("           res = "+uInFuel+" * (1 - 2*H); ");
            stw.WriteLine("      }   ");
            stw.WriteLine("  else { ");
            stw.WriteLine("      double H = 0.5 * (1.0 + Math.Tanh("+sigma+" * (X[1] + ("+inletRadius+")))); ");
            stw.WriteLine("           res = "+uInFuel+" * ( 2*H-1); ");
            stw.WriteLine("  } ");
            stw.WriteLine("return res;"); 
            stw.WriteLine("  }");


            stw.WriteLine("  static public double RegularizedPlugFlowAir(double[] X) {");         
            stw.WriteLine(" double res = 0;");
            stw.WriteLine("         if(X[1] > 0) { ");
            stw.WriteLine("     double H = 0.5 * (1.0 + Math.Tanh("+sigma+" * (X[1] - "+inletRadius+"))); ");
            stw.WriteLine("           res = "+uInAir+" * (1 - 2*H)*(-1); ");
            stw.WriteLine("      }   ");
            stw.WriteLine("  else { ");
            stw.WriteLine("      double H = 0.5 * (1.0 + Math.Tanh("+sigma+" * (X[1] + ("+inletRadius+")))); ");
            stw.WriteLine("           res = "+uInAir+" * ( 2*H-1)*(-1); ");
            stw.WriteLine("  } ");
            stw.WriteLine("return res;"); 
            stw.WriteLine("  }");


            stw.WriteLine("}"); 
            return stw.ToString();
        }
    }

    static public Formula Get_ConstantValue(double ConstVal, double inletRadius, double uInFuel, double uInAir, double sigma){
        return new Formula("BoundaryValues.ConstantValue", AdditionalPrefixCode:GetPrefixCode(ConstVal, inletRadius, uInFuel, uInAir,sigma));
    }

    static public Formula Get_ParabolaVelocityFuel(double ConstVal, double inletRadius, double uInFuel, double uInAir, double sigma){
        return new Formula("BoundaryValues.ParabolaVelocityFuel", AdditionalPrefixCode:GetPrefixCode(ConstVal, inletRadius, uInFuel, uInAir,sigma));
    }
    static public Formula Get_ParabolaVelocityAir(double ConstVal, double inletRadius, double uInFuel, double uInAir, double sigma){
        return new Formula("BoundaryValues.ParabolaVelocityAir", AdditionalPrefixCode:GetPrefixCode(ConstVal, inletRadius, uInFuel, uInAir,sigma));
    }

    static public Formula Get_RegularizedPlugFlowFuel(double ConstVal, double inletRadius, double uInFuel, double uInAir, double sigma){
        return new Formula("BoundaryValues.RegularizedPlugFlowFuel", AdditionalPrefixCode:GetPrefixCode(ConstVal, inletRadius, uInFuel, uInAir,sigma));
    }

    static public Formula Get_RegularizedPlugFlowAir(double ConstVal, double inletRadius, double uInFuel, double uInAir, double sigma){
        return new Formula("BoundaryValues.RegularizedPlugFlowAir", AdditionalPrefixCode:GetPrefixCode(ConstVal, inletRadius, uInFuel, uInAir,sigma));
    }

}

Create base control file

In this ControlFile basic configuration of the CounterDiffusionFlame is defined.

static XNSEC_Control GiveMeTheCtrlFile(int dg, int nCells, bool isMF, double massFuelIn, double massAirIn, bool parabolicVelocityProfile, bool UsefullGeometry, bool wallBounded, double xNodesMultiplier, double smoothing) {
    var CC = new ChemicalConstants();

    var C = isMF ? new XNSEC_MF_Control() : new XNSEC_Control();

    C.NumberOfChemicalSpecies = 4; 
    C.SetDGdegree(dg); //
    if(UsefullGeometry) {
        C.SetGrid(GridFactory.GenerateGrid_FullGeometry(nCells));  //
    } else {
        C.SetGrid(GridFactory.GenerateGrid(nCells, wallBounded, xNodesMultiplier));  //
    }

    C.MatParamsMode = MaterialParamsMode.Sutherland;  //
    // Problem Definition
    //===================
    double L = 0.02; // separation between the two inlets, meters
    double TemperatureInFuel = 300;  //
    double TemperatureInOxidizer = 300; //
    double AtmPressure = 101325; // Pa
    double[] FuelInletConcentrations = new double[] { 0.2, 0.0, 0.0, 0.0, 0.8 };

    double[] OxidizerInletConcentrations = new double[] { 0.0, 0.23, 0.0, 0.0, 0.77 }; //
    double[] MWs = new double[] { CC.MW_CH4, CC.MW_O2, CC.MW_CO2, CC.MW_H2O, CC.MW_N2 };

    double mwFuel = CC.getAvgMW(MWs, FuelInletConcentrations);
    double mwAir = CC.getAvgMW(MWs, OxidizerInletConcentrations);
    double densityAirIn = AtmPressure * mwAir / (CC.R_gas * TemperatureInOxidizer * 1000); // kg / m3
    double densityFuelIn = AtmPressure * mwFuel / (CC.R_gas * TemperatureInFuel * 1000); // kg / m3.
    double uInFuel = massFuelIn / densityFuelIn; //
    double uInAir = massAirIn / densityAirIn; //

    Console.WriteLine("VelocityFuel" + uInFuel);
    Console.WriteLine("VelocityAir" + uInAir);

    // Reference values
    //===================
    // Basic units to be used: Kg, m, s, mol, pa,
    double TRef = TemperatureInOxidizer;// Reference temperature  is the inlet temperature, (K)
    double pRef = AtmPressure; // Pa
    double uRef = Math.Max(uInFuel, uInAir); // m/s
    double LRef = L;
    C.GravityDirection = new double[] { 0.0, 0.0, 0.0 }; //No gravity.

    // Solver configuration
    // =======================
    C.smoothingFactor = smoothing;
    // C.NonLinearSolver.ConvergenceCriterion = 1e-8;
    // C.LinearSolver.ConvergenceCriterion = 1e-10;

    C.NonLinearSolver.verbose = true;
    C.NonLinearSolver.SolverCode = NonLinearSolverCode.Newton;
    C.NonLinearSolver.MaxSolverIterations = 10;
    C.LinearSolver = LinearSolverCode.direct_pardiso.GetConfig();
    C.TimesteppingMode = AppControl._TimesteppingMode.Steady;
    C.saveperiod = 1;
    C.PenaltyViscMomentum = 1.0; 
    C.PenaltyHeatConduction = 1.0;
    C.YFuelInlet = FuelInletConcentrations[0];
    C.YOxInlet = OxidizerInletConcentrations[1];
    C.FuelInletConcentrations = FuelInletConcentrations;
    C.OxidizerInletConcentrations = OxidizerInletConcentrations;
    C.TFuelInlet = 1.0;
    C.TOxInlet = 1.0;
    C.PhysicalParameters.IncludeConvection = true;
    // Chemical related parameters
    double s = (CC.nu_O2 * CC.MW_O2) / (CC.nu_CH4 * CC.MW_CH4);
    C.phi = s * C.YFuelInlet / C.YOxInlet;
    C.zSt = 1.0 / (1.0 + C.phi);
    var MLC = new MaterialLawCombustion(300, new double[] { }, C.MatParamsMode, C.rhoOne, true, 1.0, 1, 1, C.YOxInlet, C.YFuelInlet, C.zSt, CC, 0.75);
    var ThermoProperties = new ThermodynamicalProperties();

    //==========================
    //Derived reference values
    //==========================
    C.uRef = uRef; // Reference velocity
    C.LRef = LRef; // reference length
    C.pRef = AtmPressure; // reference pressure
    C.TRef = TemperatureInFuel;// reference temperature
    C.MWRef = MLC.getAvgMW(MWs, C.OxidizerInletConcentrations); // Air mean molecular weight
    C.rhoRef = C.pRef * C.MWRef / (8.314 * C.TRef * 1000); // Kg/m3. ok ;
    C.cpRef = 1.3;//ThermoProperties.Calculate_Cp_Mixture(new double[] { 0.23, 0.77 }, new string[] { "O2", "N2" }, 300); // 1.219185317353029;// Representative value, KJ/Kg K ========> 1.31 for the one-step kinetic model
    C.muRef = MLC.getViscosityDim(300);
    C.MolarMasses = new double[] { C.CC.MW_CH4, C.CC.MW_O2, C.CC.MW_CO2, C.CC.MW_H2O, C.CC.MW_N2 };
    C.MolarMasses.ScaleV(1.0 / C.MWRef); //NonDimensionalized Molar masses
    C.T_ref_Sutherland = 300;
    double heatRelease_Ref = (C.TRef * C.cpRef);
    C.HeatRelease = C.CC.HeatReleaseMass / heatRelease_Ref;
    C.B = CC.PreExponentialFactor;

    C.StoichiometricCoefficients = new double[] { -1, -2, 1, 2, 0 };

    C.Damk = C.rhoRef * C.LRef * C.B / (C.uRef * C.MWRef);
    C.Reynolds = C.rhoRef * C.uRef * C.LRef / C.muRef;
    C.Prandtl = 0.75; 
    C.Schmidt = C.Prandtl; // Because Lewis number  is assumed as 1.0  (Le = Pr/Sc)
    C.Lewis = new double[] { 0.97, 1.11, 1.39, 0.83, 1.0 };
    // C.Lewis = new double[] { 1.0, 1.0, 1.0, 1.0, 1.0 };

    double g = 9.8; // m/s2
    C.Froude = Math.Sqrt(uRef * uRef / (C.LRef * g)); // Not used
    C.T_ref_Sutherland = 300; //////// Check this
    C.ReactionRateConstants = new double[] { C.Damk, CC.Ta / TRef, 1.0, 1.0 }; 
    //==========================
    // Initial conditions
    //==========================
    double dummy = 0;
    double Radius = 0.5;
    C.AddInitialValue(VariableNames.VelocityX, BoundaryValueFactory.Get_ConstantValue(0.0, Radius, uInFuel / C.uRef, uInAir / C.uRef,dummy));
    C.AddInitialValue(VariableNames.VelocityY, BoundaryValueFactory.Get_ConstantValue(0.0, Radius, uInFuel / C.uRef, uInAir / C.uRef,dummy));
    C.AddInitialValue(VariableNames.Pressure, BoundaryValueFactory.Get_ConstantValue(0.0, Radius, uInFuel / C.uRef, uInAir / C.uRef,dummy));

    //==========================
    // Boundary conditions
    //==========================
    double sigma = 10;
    if(parabolicVelocityProfile) {
        C.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.Velocity_d(0), BoundaryValueFactory.Get_ParabolaVelocityFuel(dummy, Radius, uInFuel / C.uRef, dummy, dummy));
        C.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.Velocity_d(0), BoundaryValueFactory.Get_ParabolaVelocityAir(dummy, Radius, dummy, uInAir / C.uRef, dummy));
    } else { // Plug Flow
        // C.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.Velocity_d(0), BoundaryValueFactory.Get_ConstantValue(uInFuel / C.uRef, dummy, dummy, dummy, dummy));
        // C.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.Velocity_d(0), BoundaryValueFactory.Get_ConstantValue((-1) * uInAir / C.uRef, dummy, dummy, dummy, dummy));
        C.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.Velocity_d(0), BoundaryValueFactory.Get_RegularizedPlugFlowFuel(dummy, Radius, uInFuel / C.uRef, dummy, sigma));
        C.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.Velocity_d(0), BoundaryValueFactory.Get_RegularizedPlugFlowAir( dummy, Radius, dummy,  uInAir / C.uRef, sigma));
    }

    C.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.Velocity_d(1), BoundaryValueFactory.Get_ConstantValue(0.0, dummy, dummy, dummy,dummy));
    C.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.Velocity_d(1), BoundaryValueFactory.Get_ConstantValue(0.0, dummy, dummy, dummy,dummy));
    if(wallBounded){
        C.AddBoundaryValue("Wall", VariableNames.VelocityX, BoundaryValueFactory.Get_ConstantValue(0.0, dummy, dummy, dummy,dummy));
        C.AddBoundaryValue("Wall", VariableNames.VelocityY, BoundaryValueFactory.Get_ConstantValue(0.0, dummy, dummy, dummy,dummy));
    }
    return C;
}

Starting the MixtureFraction simulation

Configuration for the simulation using the mixture fraction approach, where an infinite reaction rate is assumed. Used to find adequate starting solution for the full problem.

static XNSEC_Control GiveMeTheMixtureFractionCtrlFile(int dg, int nCells, double massFuelIn, double massAirIn, bool parabolicVelocityProfile, int multiplier, bool useFullGeometry, bool wallBounded, double xNodesMultiplier, double smoothing){
    var C_MixtureFraction = GiveMeTheCtrlFile(dg, nCells, true,massFuelIn, massAirIn, parabolicVelocityProfile, useFullGeometry,  wallBounded , xNodesMultiplier, smoothing);
    C_MixtureFraction.physicsMode = PhysicsMode.MixtureFraction;
    C_MixtureFraction.ProjectName = "CounterDifFlame";
    string name = C_MixtureFraction.ProjectName + "P" + dg + "K" + nCells + "Multiplier" +multiplier + "xNodesMultiplier"+xNodesMultiplier+"smoothing"+smoothing;
    C_MixtureFraction.SessionName = "FS_" + name;

    C_MixtureFraction.UseSelfMadeTemporalOperator = false;
    C_MixtureFraction.ChemicalReactionActive = false;
    C_MixtureFraction.physicsMode = PhysicsMode.MixtureFraction;
    C_MixtureFraction.NonLinearSolver.MaxSolverIterations = 50; 
    C_MixtureFraction.dummycounter = multiplier;
    
    // Boundary and initial conditions
    double dummy = -11111111;
    C_MixtureFraction.AddInitialValue(VariableNames.MixtureFraction,BoundaryValueFactory.Get_ConstantValue(1.0,dummy,dummy , dummy, dummy));
    C_MixtureFraction.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.MixtureFraction, BoundaryValueFactory.Get_ConstantValue(1.0,dummy,dummy , dummy, dummy));
    C_MixtureFraction.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.MixtureFraction, BoundaryValueFactory.Get_ConstantValue(0.0,dummy,dummy , dummy, dummy));

    C_MixtureFraction.Tags.Add(multiplier.ToString()); // Tag used for restart of the full problem
    double radius_inlet = 0.5; // radius is 1/2 from length separation
    var troubledPoints = new List<double[]>();
    if(useFullGeometry){
        troubledPoints.Add(new double[]{2, +radius_inlet });
        troubledPoints.Add(new double[]{2, -radius_inlet });
        troubledPoints.Add(new double[]{1, +radius_inlet });
        troubledPoints.Add(new double[]{1, -radius_inlet });
    } else {
        troubledPoints.Add(new double[]{0, +radius_inlet });
        troubledPoints.Add(new double[]{0, -radius_inlet });
        troubledPoints.Add(new double[]{1, +radius_inlet });
        troubledPoints.Add(new double[]{1, -radius_inlet });
    }
    

    bool useHomotopy = false;

    if(useHomotopy) { 
        C_MixtureFraction.HomotopyApproach = XNSEC_Control.HomotopyType.Automatic;
        // C_MixtureFraction.HomotopyVariable = XNSEC_Control.HomotopyVariableEnum.VelocityInletMultiplier;
        // C_MixtureFraction.homotopieAimedValue = multiplier;
        C_MixtureFraction.HomotopyVariable = XNSEC_Control.HomotopyVariableEnum.Reynolds;
        C_MixtureFraction.homotopieAimedValue = C_MixtureFraction.Reynolds;
    }


    C_MixtureFraction.AdaptiveMeshRefinement = true;
    C_MixtureFraction.AMR_startUpSweeps =2;
    int NoOfPseudoTimesteps = 2;
    C_MixtureFraction.TimesteppingMode = BoSSS.Solution.Control.AppControl._TimesteppingMode.Steady; 
    C_MixtureFraction.NoOfTimesteps = NoOfPseudoTimesteps ;

    // C_MixtureFraction.activeAMRlevelIndicators.Add( new BoSSS.Application.XNSEC.AMR_onProblematicPoints(troubledPoints,C_MixtureFraction.AMR_startUpSweeps) );  
    C_MixtureFraction.activeAMRlevelIndicators.Add(new BoSSS.Application.XNSEC.AMR_RefineAroundProblematicPoints(troubledPoints, 1, 0.01));
    C_MixtureFraction.activeAMRlevelIndicators.Add( new BoSSS.Application.XNSEC.AMR_onFlameSheet(C_MixtureFraction.zSt,3) );
    return C_MixtureFraction;
}

Send and run jobs

foreach(double xDensityMult in  xCellsMultiplier)
foreach(int nCells in nCellsArray){

foreach(int i in multiplierS){
  foreach(double smooth in Smoothings) {
  double massFuelIn = InitialMassFuelIn*i;
  double massAirIn = InitialMassAirIn*i;

  Type solver_MF = typeof(BoSSS.Application.XNSEC.XNSEC_MixtureFraction);
  var C_MixtureFraction = GiveMeTheMixtureFractionCtrlFile(dgMF, nCells, massFuelIn, massAirIn,  parabolicVelocityProfile, i, useFullGeometry, wallBounded, xDensityMult,smooth);
  string jobName       = C_MixtureFraction.SessionName;
  Console.WriteLine(jobName);
  var oneJob           = new Job(jobName, solver_MF);
  oneJob.NumberOfMPIProcs =  8;
  oneJob.SetControlObject(C_MixtureFraction);
    // oneJob.Activate(myBatch); 
    oneJob.Activate(); 
}
}
}

Grid Edge Tags changed.

An equivalent grid (e874f3f1-fae4-4a7e-b30e-6096257af25e) is already present in the database -- the grid will not be saved.

VelocityFuel0.04852837051919214

VelocityAir0.12295612139156184

FS_CounterDifFlameP3K10Multiplier2xNodesMultiplier3smoothing80

Deployments so far (1): (Job token: 802872, FinishedSuccessful 'CounterFlowFlame_MF_FullComparison-XNSEC2023Oct30_104008.531442' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "CounterFlowFlame_MF_FullComparison	FS_CounterDifFlameP3K10Multiplier2xNodesMultiplier3smoothing80	10/30/2023 10:40:30	d35146c1..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Grid Edge Tags changed.

An equivalent grid (e874f3f1-fae4-4a7e-b30e-6096257af25e) is already present in the database -- the grid will not be saved.

VelocityFuel0.12132092629798034

VelocityAir0.3073903034789046

FS_CounterDifFlameP3K10Multiplier5xNodesMultiplier3smoothing80

Deployments so far (1): (Job token: 802874, FinishedSuccessful 'CounterFlowFlame_MF_FullComparison-XNSEC2023Oct30_104026.269782' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "CounterFlowFlame_MF_FullComparison	FS_CounterDifFlameP3K10Multiplier5xNodesMultiplier3smoothing80	10/30/2023 10:40:49	bb35031d..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Grid Edge Tags changed.

An equivalent grid (e874f3f1-fae4-4a7e-b30e-6096257af25e) is already present in the database -- the grid will not be saved.

VelocityFuel0.2669060378555568

VelocityAir0.6762586676535901

FS_CounterDifFlameP3K10Multiplier11xNodesMultiplier3smoothing80

Deployments so far (1): (Job token: 802877, FinishedSuccessful 'CounterFlowFlame_MF_FullComparison-XNSEC2023Oct30_104045.116683' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "CounterFlowFlame_MF_FullComparison	FS_CounterDifFlameP3K10Multiplier11xNodesMultiplier3smoothing80	10/30/2023 10:41:10	09e6449e..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

BoSSSshell.WorkflowMgm.BlockUntilAllJobsTerminate();

All jobs finished.

Starting the finite-rate chemistry simulation

Now that the simulation for an "infinite" reaction rate is done, we use it for initializing the system with finite reaction rate. The goal is to obtain solutions of the counter difussion flame for increasing strain values. We start with a low strain (bigger Dahmkoehler number), which is increased until extintion is (hopefully) found

static XNSEC_Control GiveMeTheFullCtrlFile(int dg, int nCells, double massFuelIn, double massAirIn, ISessionInfo SessionToRestart, bool parabolicVelocityProfile, bool chemReactionActive, bool useFullGeometry, bool wallBounded, int mult, int counter, double xNodesMultiplier, double smoothing) {
    var C_OneStep = GiveMeTheCtrlFile(dg, nCells, false, massFuelIn, massAirIn, parabolicVelocityProfile, useFullGeometry, wallBounded, xNodesMultiplier, smoothing);
    C_OneStep.physicsMode = PhysicsMode.Combustion;
    C_OneStep.ProjectName = "CounterDifFlame";
    string name = C_OneStep.ProjectName + "P" + dg + "K" + nCells + "mult" + mult+"_c"+counter ;
    C_OneStep.SessionName = "Full_" + name;
    C_OneStep.VariableOneStepParameters = true;

    // C_OneStep.Tags.Add("VelocityMultiplier" + mult);
    C_OneStep.Tags.Add( mult.ToString());
    C_OneStep.dummycounter = counter;
    C_OneStep.UseSelfMadeTemporalOperator = false;
    C_OneStep.myThermalWallType = SIPDiffusionTemperature.ThermalWallType.Adiabatic;
    C_OneStep.Timestepper_LevelSetHandling = BoSSS.Solution.XdgTimestepping.LevelSetHandling.None;
    C_OneStep.UseMixtureFractionsForCombustionInitialization = true;

    C_OneStep.LinearSolver = new BoSSS.Solution.AdvancedSolvers.OrthoMGSchwarzConfig() {
        NoOfMultigridLevels = 5,
        //verbose = true
    };

    C_OneStep.ChemicalReactionActive = chemReactionActive;
    C_OneStep.AdaptiveMeshRefinement = true;
    
    C_OneStep.HeatCapacityMode = MaterialLaw_MultipleSpecies.CpCalculationMode.mixture;

    bool AMRinEachNewtonStep = false;

    if( AMRinEachNewtonStep) {
        C_OneStep.NoOfTimesteps = 4;
        C_OneStep.NonLinearSolver.MaxSolverIterations = 8; // Do only one newton iteration before refining
        C_OneStep.NonLinearSolver.MinSolverIterations = 8; // Do only one newton iteration before refining
    } else{
        C_OneStep.NoOfTimesteps = 3; // The steady solution will be calculated again and do AMR
        C_OneStep.NonLinearSolver.MaxSolverIterations = 3;
    }

    C_OneStep.AMR_startUpSweeps = 0;
    if(C_OneStep.ChemicalReactionActive){
        C_OneStep.activeAMRlevelIndicators.Add(new AMR_onReactiveZones( 2, 0.2));
        C_OneStep.activeAMRlevelIndicators.Add(new AMR_BasedOnVariableLimits("Temperature", new double[] { -100, 4 },3)); // Refine all cells with T > 5 (and T < -100)
        // C_OneStep.activeAMRlevelIndicators.Add(new AMR_BasedOnFieldGradient(2, 0.9, VariableNames.Temperature));
        // C_OneStep.activeAMRlevelIndicators.Add(new AMR_BasedOnPerssonSensor(VariableNames.Temperature, 2));
    }
    // C_OneStep.activeAMRlevelIndicators.Add(new AMR_BasedOnPerssonSensor(VariableNames.Temperature, 3));
    // C_OneStep.NonLinearSolver.MaxSolverIterations = 10;

    // limiting of variable values
    Dictionary<string, Tuple<double, double>> Bounds = new Dictionary<string, Tuple<double, double>>();
    double eps = 1e-2;
    Bounds.Add(VariableNames.Temperature, new Tuple<double, double>(1.0 - eps, 10)); // Min temp should be the inlet temperature.
    Bounds.Add(VariableNames.MassFraction0, new Tuple<double, double>(0.0 - 1e-4, 1.0 + 1e-4)); // Between 0 and 1 per definition
    Bounds.Add(VariableNames.MassFraction1, new Tuple<double, double>(0.0 - 1e-4, 1.0 + 1e-4));
    Bounds.Add(VariableNames.MassFraction2, new Tuple<double, double>(0.0 - 1e-4, 1.0 + 1e-4));
    Bounds.Add(VariableNames.MassFraction3, new Tuple<double, double>(0.0 - 1e-4, 1.0 + 1e-4));
    C_OneStep.VariableBounds = Bounds;
    // Boundary conditions

    double dummy = 0;
    if(SessionToRestart != null) {
        C_OneStep.SetRestart(SessionToRestart);
    } else {
        C_OneStep.AddInitialValue(VariableNames.Temperature, BoundaryValueFactory.Get_ConstantValue(1.0, dummy, dummy, dummy, dummy));
        C_OneStep.AddInitialValue(VariableNames.MassFraction0, BoundaryValueFactory.Get_ConstantValue(0.0, dummy, dummy, dummy, dummy));
        C_OneStep.AddInitialValue(VariableNames.MassFraction1, BoundaryValueFactory.Get_ConstantValue(0.23, dummy, dummy, dummy, dummy));
        C_OneStep.AddInitialValue(VariableNames.MassFraction2, BoundaryValueFactory.Get_ConstantValue(0.0, dummy, dummy, dummy, dummy));
        C_OneStep.AddInitialValue(VariableNames.MassFraction3, BoundaryValueFactory.Get_ConstantValue(0.0, dummy, dummy, dummy, dummy));
    }
    if(wallBounded){
        C_OneStep.AddBoundaryValue("Wall", VariableNames.Temperature, BoundaryValueFactory.Get_ConstantValue(1.0, dummy, dummy, dummy, dummy));
    }
    C_OneStep.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.Temperature, BoundaryValueFactory.Get_ConstantValue(1.0, dummy, dummy, dummy, dummy));
    C_OneStep.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.MassFraction0, BoundaryValueFactory.Get_ConstantValue(C_OneStep.FuelInletConcentrations[0], dummy, dummy, dummy, dummy));
    C_OneStep.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.MassFraction1, BoundaryValueFactory.Get_ConstantValue(C_OneStep.FuelInletConcentrations[1], dummy, dummy, dummy, dummy));
    C_OneStep.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.MassFraction2, BoundaryValueFactory.Get_ConstantValue(C_OneStep.FuelInletConcentrations[2], dummy, dummy, dummy, dummy));
    C_OneStep.AddBoundaryValue("Velocity_Inlet_CH4", VariableNames.MassFraction3, BoundaryValueFactory.Get_ConstantValue(C_OneStep.FuelInletConcentrations[3], dummy, dummy, dummy, dummy));

    C_OneStep.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.Temperature, BoundaryValueFactory.Get_ConstantValue(1.0, dummy, dummy, dummy, dummy));
    C_OneStep.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.MassFraction0, BoundaryValueFactory.Get_ConstantValue(C_OneStep.OxidizerInletConcentrations[0], dummy, dummy, dummy, dummy));
    C_OneStep.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.MassFraction1, BoundaryValueFactory.Get_ConstantValue(C_OneStep.OxidizerInletConcentrations[1], dummy, dummy, dummy, dummy));
    C_OneStep.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.MassFraction2, BoundaryValueFactory.Get_ConstantValue(C_OneStep.OxidizerInletConcentrations[2], dummy, dummy, dummy, dummy));
    C_OneStep.AddBoundaryValue("Velocity_Inlet_O2", VariableNames.MassFraction3, BoundaryValueFactory.Get_ConstantValue(C_OneStep.OxidizerInletConcentrations[3], dummy, dummy, dummy, dummy));
    return C_OneStep;
}

Type solver = typeof(BoSSS.Application.XNSEC.XNSEC);

Calculate the full solution for the initial value

int counter = 0;
foreach (double xDensityNodes in xCellsMultiplier){
    foreach (int nCells in nCellsArray) {
        foreach (int dg in DGDegrees) {
            foreach(double smooth in Smoothings){
            foreach (int i in multiplierS) {
                var sess = (myDb.Sessions.Where(s => s.Name == "FS_CounterDifFlameP" + dgMF + "K" + nCells + "Multiplier" + i + "xNodesMultiplier" + xDensityNodes+"smoothing"+smooth)).FirstOrDefault();
                var C = GiveMeTheFullCtrlFile(dg, nCells, InitialMassFuelIn * i, InitialMassAirIn * i, sess, parabolicVelocityProfile, chemicalReactionActive, useFullGeometry, wallBounded, i, counter, xDensityNodes,smooth);
                string jobName = C.SessionName;
                Console.WriteLine(jobName);
                var oneJob = new Job(jobName, solver);
                oneJob.NumberOfMPIProcs = 12;
                oneJob.SetControlObject(C);
                  // oneJob.Activate(myBatch); 
                oneJob.Activate(); 
                counter++;
            }
        }
    }
}
}

Grid Edge Tags changed.

An equivalent grid (e874f3f1-fae4-4a7e-b30e-6096257af25e) is already present in the database -- the grid will not be saved.

VelocityFuel0.04852837051919214

VelocityAir0.12295612139156184

Full_CounterDifFlameP4K10mult2_c0

Deployments so far (1): (Job token: 802973, FinishedSuccessful 'CounterFlowFlame_MF_FullComparison-XNSEC2023Oct30_110004.794238' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "CounterFlowFlame_MF_FullComparison	Full_CounterDifFlameP4K10mult2_c0	10/30/2023 11:00:29	43b71a34..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Grid Edge Tags changed.

An equivalent grid (e874f3f1-fae4-4a7e-b30e-6096257af25e) is already present in the database -- the grid will not be saved.

VelocityFuel0.12132092629798034

VelocityAir0.3073903034789046

Full_CounterDifFlameP4K10mult5_c1

Deployments so far (1): (Job token: 802974, FinishedSuccessful 'CounterFlowFlame_MF_FullComparison-XNSEC2023Oct30_110025.790218' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "CounterFlowFlame_MF_FullComparison	Full_CounterDifFlameP4K10mult5_c1	10/30/2023 11:00:51	3638c5e8..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

Grid Edge Tags changed.

An equivalent grid (e874f3f1-fae4-4a7e-b30e-6096257af25e) is already present in the database -- the grid will not be saved.

VelocityFuel0.2669060378555568

VelocityAir0.6762586676535901

Full_CounterDifFlameP4K10mult11_c2

Deployments so far (1): (Job token: 802976, FinishedSuccessful 'CounterFlowFlame_MF_FullComparison-XNSEC2023Oct30_110046.431066' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "CounterFlowFlame_MF_FullComparison	Full_CounterDifFlameP4K10mult11_c2	10/30/2023 11:01:12	37f41a82..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

// wait for all jobs to finish (up to 2 days, check every 10 minutes)
BoSSSshell.WorkflowMgm.BlockUntilAllJobsTerminate(TimeOutSeconds:(3600*24*2), PollingIntervallSeconds:(60*10));

All jobs finished.

// detect failed Jobs in the job management
var suspects = BoSSSshell.WorkflowMgm.AllJobs.Select(kv => kv.Value)
    .Where(job => job.LatestSession.Tags.Contains(SessionInfo.NOT_TERMINATED_TAG)
                  || job.LatestSession.Tags.Contains(SessionInfo.SOLVER_ERROR)).ToArray();
suspects

NUnit.Framework.Assert.IsTrue(suspects.Count() <= 0, $"{suspects.Count()} Failed Jobs of {BoSSSshell.WorkflowMgm.AllJobs.Count()} in total.");

// // Delete all MF calculations
// Console.WriteLine("Deleting mass fraction calculations");
// myDb.Sessions.Where(sess => sess.Name.Contains("FS")).ForEach(sess=> sess.Delete(true))