// #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.

int Resolution = 9; // was 9 before for all calculations
int dg = 2;
double[] ExpansionRatioS = new double[]{2.0};
bool[] varParams = new bool[]{false};

BoSSSshell.WorkflowMgm.Init("HeatedBackwardFacingStep"); // HeatedBackwardFacingStep

Project name is set to 'HeatedBackwardFacingStep'.

Default Execution queue is chosen for the database.

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

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

double Sdim = 15;  //mm, dim step height
double hdim = Sdim/(2-1.0);
double Hdim = hdim+Sdim; // mm channel height


double LREF = Sdim;// step height, paper from XIE

// static double LREF = 2*hdim; // Hydraulic referemce, paper from Armaly
double H = Hdim / LREF;// nondim channel height
double S = Sdim / LREF;  //nondim step height
double h = hdim/LREF;
double L0 = S; // nondim upstream wall
double L = 50 * S; // downstream bc

Create grid

public static class GridFactory {

    public static Grid2D GenerateGrid(int Res, double S, double h, double H, double L0,double L) {
        // Grid2D grd;

        double[] CutOut1Point1 = new double[2] { 0.0, 0.0 };
        double[] CutOut1Point2 = new double[2] { -L0, S };

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

        
        // left part
        var _leftNodes = GenericBlas.SinLinSpacing(-2 * L0, 0, 0.99, (2 * Res) + 1);
        var leftNodes = _leftNodes.GetSubVector(_leftNodes.Length / 2, _leftNodes.Length / 2 + 1);
        //right part
        // int xNodesMultiplier  = 110;
        int xNodesMultiplier  = 230;
        // int xNodesMultiplier  = 50;
        // int xNodesMultiplier  = 5;

        var _rightNodes = GenericBlas.SinLinSpacing(0, 2 * L, 0.9, (xNodesMultiplier * Res) + 1); 
        var rightNodes = _rightNodes.GetSubVector(0, _rightNodes.Length / 2 + 1);

        List<double> listX = new List<double>();
        listX.AddRange(leftNodes.Take(leftNodes.Count() - 1).ToList());
        listX.AddRange(rightNodes.Take(rightNodes.Count() + 1).ToList());
        var xNodes = listX.ToArray();


        //Bottom part
        var _bottomNodes = GenericBlas.SinLinSpacing(0, S, 0.8, 2*(2 * Res) + 1); 
        //upper part
        var _upperNodes = GenericBlas.SinLinSpacing(S, S+h, 0.8, 2*(2 * Res) + 1);
        // var upperNodes = _upperNodes.GetSubVector(0, _upperNodes.Length / 2 + 1);        
        // for(int i = 0; i < upperNodes.Count(); i++){
        //     upperNodes[i] = upperNodes[i] + S;
        // }
        
        List<double> listY = new List<double>();
        listY.AddRange(_bottomNodes.Take(_bottomNodes.Count() - 1).ToList());
        listY.AddRange(_upperNodes.Take(_upperNodes.Count() + 1).ToList());
        var yNodes = listY.ToArray();
        var grd = Grid2D.Cartesian2DGrid(xNodes, yNodes, CutOuts: CutOut1);

        Console.WriteLine("Number of nodes in X direction:" +listX.Count() );
        Console.WriteLine("Number of nodes in Y direction:" +yNodes.Count() );
        grd.EdgeTagNames.Add(1, "Velocity_Inlet");
        grd.EdgeTagNames.Add(3, "Pressure_Outlet");
        grd.EdgeTagNames.Add(2, "NoSlipNeumann");
        grd.EdgeTagNames.Add(4, "Wall");

        grd.DefineEdgeTags(delegate (double[] X) {
            double x = X[0];
            double y = X[1];

            //Inlet oxidizer
            if (Math.Abs(x + L0) < 1e-8)
                return 1;

            //Outlet
            if (Math.Abs(x - L) < 1e-8)
                return 3;

            //Bottom Wall
            if (Math.Abs(y - 0) < 1e-8 && x > 1e-8)
                return 4;

            return 2; 
        });

        myDb.SaveGrid(ref grd);

        return grd;
    }
}

public static class BoundaryValueFactory { 

    public static string GetPrefixCode(double const_val, double S, double h, double umean ) {
        using(var stw = new System.IO.StringWriter()) {
           
           stw.WriteLine("static class BoundaryValues {");


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

           double H = S+h;           
           stw.WriteLine("  static public double ParabolaVelocity(double[] X) {");
           stw.WriteLine("    return  -6 * "+umean+" * (X[1] -"+ S+") * (X[1] - 1 * "+H+") / ("+h * h+") ;");
           stw.WriteLine("  }");


           stw.WriteLine("}"); 
           return stw.ToString();
        }
    }
    static public Formula Get_ConstantVal(double const_val, double S, double h,double umean){
        return new Formula("BoundaryValues.ConstantVal", AdditionalPrefixCode:GetPrefixCode(const_val,S,h,umean));
    }

    static public Formula Get_ParabolaVelocity(double const_val, double S, double h, double umean){
        return new Formula("BoundaryValues.ParabolaVelocity", AdditionalPrefixCode:GetPrefixCode(const_val,S,h,umean));
    }

}

Send and run jobs

static XNSEC_Control GiveMeTheCtrlFile(int dg, int Resolution, bool varParams, double expansionRatio) {
    XNSEC_Control C = new XNSEC_Control();
    Console.WriteLine("////////////////////////////////////////////////////////////////////////////////////");
    Console.WriteLine("BackwardFacingStep2");
    Console.WriteLine("////////////////////////////////////////////////////////////////////////////////////");

    // Solver configuration
    // ==============

    C.TimesteppingMode = AppControl._TimesteppingMode.Steady;

    C.NonLinearSolver.SolverCode = NonLinearSolverCode.Newton;
    C.NonLinearSolver.verbose = true;
    C.ProjectName = "Heated_BackwardFacingStep";
    C.SessionName = C.ProjectName + "dg" + dg + "k" + Resolution +  "ER" + expansionRatio+ "VP"+varParams.ToString() ;
    //C.physicsMode = PhysicsMode.LowMach;
    C.physicsMode = PhysicsMode.Combustion;
    C.EnableMassFractions = false;
    C.EnableTemperature = true;
    
    C.Tags.Add(expansionRatio.ToString()); //
    C.Tags.Add(varParams.ToString());

   


    if (varParams) {
        C.MatParamsMode = MaterialParamsMode.Sutherland;
        C.rhoOne = false;
    } else{
        C.rhoOne = true;
        C.MatParamsMode = MaterialParamsMode.Constant;
    }
    C.AnalyticsolutionSwitch = false;

    C.PhysicalParameters.IncludeConvection = true;

    C.NumberOfChemicalSpecies = 1;
    C.ChemicalReactionActive = false;
    C.SetDGdegree(dg);
    C.GravityDirection = new double[] { 0.0, 0.0, 0.0 };
    C.HomotopyApproach = XNSEC_Control.HomotopyType.Manual;
    C.HomotopyVariable = XNSEC_Control.HomotopyVariableEnum.Reynolds;  

    C.HomotopyArray = new double[]{1e-2,1e-1,1e0,10,50,100,150,200,250,300,350,400,450,500,550,600,650,700};

    C.HeatRelease = 0.0;
    // Parameters
    // ==============

    C.Reynolds = 700;
    C.Prandtl = 0.7132;
    C.Schmidt = 1.0;
    C.PenaltyViscMomentum = 1.0;
    C.PenaltyHeatConduction = 1.0;

    C.NonLinearSolver.verbose = true;
    
    // C.NonLinearSolver.ConvergenceCriterion = 1e-10;
    C.NonLinearSolver.MaxSolverIterations = 30;
    C.LinearSolver = LinearSolverCode.direct_pardiso.GetConfig();
    // C.LinearSolver.NoOfMultigridLevels = 7;
    // C.TimesteppingMode = AppControl._TimesteppingMode.Transient;
    // C.NoOfTimesteps = 10000;
    // C.dtFixed = 0.5;
    // C.TracingNamespaces = "BoSSS.Solution.AdvancedSolvers.Newton"; ////////////////////////////////////////////////////////////////////////////////////////////////////////

    // C.HomotopyApproach = XNSEC_Control.HomotopyType.Automatic;
    // C.HomotopyVariable = XNSEC_Control.HomotopyVariableEnum.Reynolds;
    // C.homotopieAimedValue = C.Reynolds;
    // C.StartingHomotopyValue = 100; // Suficiently easy to find solution
    // C.NonLinearSolver.HomotopyStepLongFail = 10;

    // Grid declaration
    // ===============

     double Sdim = 15;  //mm, dim step height
     double hdim = Sdim/(expansionRatio-1.0);
     double Hdim = hdim+Sdim; // mm channel height


     double LREF = Sdim;// step height, paper from XIE

    // static double LREF = 2*hdim; // Hydraulic referemce, paper from Armaly
     double H = Hdim / LREF;// nondim channel height
     double S = Sdim / LREF;  //nondim step height
     double h = hdim/LREF;
     double L0 = S; // nondim upstream wall
    double L = 50 * S; // downstream bc

    

    C.SetGrid(GridFactory.GenerateGrid(Resolution, S,h,H,L0,L));  //

    double TRef = 273 + 10;
    double Tin = (273 + 10) / TRef;
    double Twall = (273 + 40) / TRef;

    // initial values
    // ==============

    double dummy = 0;
    C.AddInitialValue(VariableNames.VelocityX, BoundaryValueFactory.Get_ConstantVal(0.5, dummy, dummy, dummy));
    C.AddInitialValue(VariableNames.VelocityY, BoundaryValueFactory.Get_ConstantVal(0.0, dummy, dummy, dummy));
    C.AddInitialValue(VariableNames.Pressure, BoundaryValueFactory.Get_ConstantVal(0.0, dummy, dummy, dummy));
    C.AddInitialValue(VariableNames.Temperature, BoundaryValueFactory.Get_ConstantVal((Tin+Twall)*0.5, dummy, dummy, dummy));
    C.AddInitialValue(VariableNames.MassFraction0, BoundaryValueFactory.Get_ConstantVal(1.0, dummy, dummy, dummy));

    // boundary conditions
    // ===================

    C.AddBoundaryValue("Wall", VariableNames.Temperature + "#A", BoundaryValueFactory.Get_ConstantVal(Twall, dummy, dummy, dummy));
    C.AddBoundaryValue("Wall", VariableNames.MassFraction0 + "#A", BoundaryValueFactory.Get_ConstantVal(1.0, dummy, dummy, dummy));



    double velMean= 1.0;
    C.AddBoundaryValue("Velocity_Inlet", VariableNames.Velocity_d(0) + "#A", BoundaryValueFactory.Get_ParabolaVelocity(-11111.0, S,h, velMean));
    C.AddBoundaryValue("Velocity_Inlet", VariableNames.Velocity_d(1) + "#A", BoundaryValueFactory.Get_ConstantVal(0.0, dummy, dummy, dummy));
    C.AddBoundaryValue("Velocity_Inlet", VariableNames.Temperature + "#A", BoundaryValueFactory.Get_ConstantVal(Tin, dummy, dummy, dummy));
    C.AddBoundaryValue("Velocity_Inlet", VariableNames.MassFraction0 + "#A", BoundaryValueFactory.Get_ConstantVal(1.0, dummy, dummy, dummy));

    C.AddBoundaryValue("Pressure_Outlet");
    return C;
}

Run Simulations

foreach(double ER in ExpansionRatioS){
        foreach(bool b in varParams){
        Type solver = typeof(BoSSS.Application.XNSEC.XNSEC);

        var C = GiveMeTheCtrlFile(dg, Resolution,  b, ER); 
        string jobName       = C.SessionName;
        Console.WriteLine(jobName);
        var oneJob           = new Job(jobName, solver);
        oneJob.SetControlObject(C);
        oneJob.NumberOfMPIProcs = MPI_PROCS;
        oneJob.Activate(myBatch); 
    }
}

////////////////////////////////////////////////////////////////////////////////////

BackwardFacingStep2

////////////////////////////////////////////////////////////////////////////////////





Number of nodes in X direction:1045

Number of nodes in Y direction:73

Grid Edge Tags changed.

An equivalent grid (c08618c1-2ea2-4d8d-97f2-b9af9fc6e7cc) is already present in the database -- the grid will not be saved.

Heated_BackwardFacingStepdg2k9ER2VPFalse

Deployments so far (1): (Job token: 802879, FinishedSuccessful 'HeatedBackwardFacingStep-XNSEC2023Oct30_104059.163837' @ MS HPC client  MSHPC-AllNodes @DC2, @\\fdygitrunner\ValidationTests\deploy, FinishedSuccessful);

Success: 1

Info: Found successful session "HeatedBackwardFacingStep	Heated_BackwardFacingStepdg2k9ER2VPFalse	10/30/2023 10:41:33	d997fb70..." -- job is marked as successful, no further action.

No submission, because job status is: FinishedSuccessful

BoSSSshell.WorkflowMgm.BlockUntilAllJobsTerminate();

All jobs finished.

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

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.");