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

using BoSSS.Foundation.Quadrature;
using BoSSS.Foundation.XDG;
using MathNet.Numerics;
using MathNet.Numerics.Interpolation;
using MathNet.Numerics.RootFinding;

BoSSSshell.WorkflowMgm.Init("HeatedBackwardFacingStep");

Project name is set to 'HeatedBackwardFacingStep'.

Default Execution queue is chosen for the database.

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

var myDb = BoSSS.Application.BoSSSpad.BoSSSshell.WorkflowMgm.DefaultDatabase;

Creation of a spline along desired edge. Used later for interpolation of variables.

// Creates a Spline on direction "d" at a given edge "em" of the field "field"
static public CubicSpline SplineOnEdge(EdgeMask em, DGField field, out double lower_Bound, out double upper_Bound, double offset = 0.0, int d = 0){

    var grd = field.GridDat;

    List<double> nodes = new List<double>();
    List<double> values = new List<double>();

    EdgeQuadrature.GetQuadrature(new int[] { 1 }, grd,
        // (new EdgeQuadratureScheme(true, em)).Compile(grd, 0),
        (new EdgeQuadratureScheme(true, em)).Compile(grd, field.Basis.Degree * 2),
        delegate (int i0, int Length, QuadRule QR, MultidimensionalArray EvalResult) {        

            MultidimensionalArray DummyIN = MultidimensionalArray.Create(Length, QR.NoOfNodes);
            MultidimensionalArray DummyOT = MultidimensionalArray.Create(Length, QR.NoOfNodes);        

            MultidimensionalArray GlobalNodes = MultidimensionalArray.Create(Length, QR.NoOfNodes, 2);        
            
            if(field is XDGField xField){
                xField.GetSpeciesShadowField("A").EvaluateEdge(i0, Length, QR.Nodes, DummyIN, DummyOT, null, null, null, null, 0, 0.0); 
            } else{
                field.EvaluateEdge(i0, Length, QR.Nodes, DummyIN, DummyOT, null, null, null, null, 0, 0.0); 
            }  

            for(int i = 0; i < Length; i++){  
                int iTrafo = ((GridData)grd).Edges.Edge2CellTrafoIndex[i0+i, 0];
                NodeSet volNodeSet = QR.Nodes.GetVolumeNodeSet(grd, iTrafo, false);
                int jCell = ((GridData)grd).Edges.CellIndices[i0+i, 0];
                grd.TransformLocal2Global(volNodeSet, jCell, 1, GlobalNodes, i);          
                int K = QR.NoOfNodes;
                for(int k = 0; k < K; k++){
                    nodes.Add(GlobalNodes[i, k, d]);
                    values.Add(DummyIN[i, k] - offset);
                }
            }

        },delegate (int i0, int Length, MultidimensionalArray ResultsOfIntegration) {
        }).Execute(); 

    lower_Bound = nodes.Min();
    upper_Bound = nodes.Max();

    return CubicSpline.InterpolateAkima(nodes.ToArray(), values.ToArray());
}

Method for finding the position where a specified field assumes a given value using the bisection method.

static public double[] PosOfValueOnEdge(double value, EdgeMask em, DGField field, int d = 0) {
    var Spline = SplineOnEdge(em, field, out double lB, out double uB, value, d);
    int numberofintervals = 30;
    double[] intervalnodes = GenericBlas.Linspace(lB, uB, numberofintervals);
    List<double> roots = new List<double>();
    for (int i = 0; i < numberofintervals - 1; i++) {
        double root = -1111;
        bool found = Bisection.TryFindRoot(t => Spline.Interpolate(t), intervalnodes[i], intervalnodes[i + 1], 1e-12, 1000, out root);
        if (found)
            roots.Add(root);
    }
    return roots.ToArray();
}

Calculation of Local Nusselt Number $\left(\text{Nu} := \frac{S}{T0-T1}\vec{n}\cdot \nabla T \right)$ along a given surface.

static public double[] LocalNusseltNumber(ITimestepInfo ts, double[] coordinates, double deltaT, int dirGradient, int dirSpline) {
    DGField Temperature       = ts.Fields.Where(f => f.Identification == "Temperature").SingleOrDefault();
    var grd = (GridData)Temperature.GridDat;
    var GradT_d = (DGField)Temperature.Clone();
    GradT_d.Clear();
    GradT_d.Derivative(1.0, Temperature, dirGradient);

    if(dirGradient == dirSpline){
    throw new Exception();
}
    EdgeMask em = new EdgeMask(grd, X => X[1] + 0.0 < 1e-12); // edge hardcoded
    var spline = SplineOnEdge(em, GradT_d, out double _a, out double  _b, d:dirSpline);
    
    double[] Nusselt = new double[coordinates.Length];
    for(int i = 0; i < coordinates.Length; i++){
        double GradT = spline.Interpolate(coordinates[i]);    
        Nusselt[i] = -1.0*GradT/deltaT;
    }   
    return Nusselt;
}

Calculation of Local Darcy frictionfactor $\text{f}_D := -\frac{8\nu}{(\hat{u}/2)^2}\vec{n}\cdot \nabla u_1$

static public double[] LocalDarcyFactor(ITimestepInfo ts, double[] coordinates, int dirGradient, int dirSpline, double Re ) {
    DGField VelocityX       = ts.Fields.Where(f => f.Identification == "VelocityX").SingleOrDefault();

        
    double nu = 15.52e-6;// Kinematic viscosity of air at 25°c, m2/s
    double S = 15.0/1000; // Reference length according to the paper of XIA, meters
    double vel = Re*nu/S; // mean velocity, m/s

    //Calculation of the velocity gradient normal to the surface
    var grd = (GridData)VelocityX.GridDat;
    var GradVelX_d = (DGField)VelocityX.Clone();
    GradVelX_d.Clear(); 
    GradVelX_d.Derivative(1.0 , VelocityX, dirGradient);

    if(dirGradient == dirSpline){
    throw new Exception();
}
    EdgeMask em = new EdgeMask(grd, X => X[1] + 0.0 < 1e-12); //  bottom wall
    var spline = SplineOnEdge(em, GradVelX_d, out double _a, out double  _b, d:dirSpline);

    double factor = 4*nu/(vel*vel/4.0);

    double[] fd = new double[coordinates.Length];
    for(int i = 0; i < coordinates.Length; i++){
        double GradVelX_dVal = spline.Interpolate(coordinates[i]);    
        fd[i] = factor*GradVelX_dVal*100 ;
    }   
    return fd;
}

This method is used for finding the position(s) of deatachment at a wall. Deattachment (or reattachment) occurs exactly at the point(s) where the friction factor is equal to zero $f_d(\vec{x})=0$

static public double[] GetPositionOfReattachment(ITimestepInfo ts, int dirGradient, int dirSpline, double ER, bool top ) {
    DGField VelocityX       = ts.Fields.Where(f => f.Identification == "VelocityX").SingleOrDefault();
    var grd = (GridData)VelocityX.GridDat;
    var GradVelX_d = (DGField)VelocityX.Clone();
    GradVelX_d.Clear();
    GradVelX_d.Derivative(1.0, VelocityX, dirGradient);

    double S = 1.0;
    double H = S+S/(ER-1.0);
    double ywall = top == true? H : 0.0;
    EdgeMask em = new EdgeMask(grd, X => Math.Abs(X[1] - ywall) < 1e-12); 
    double[] x =   PosOfValueOnEdge(0.0, em, GradVelX_d, d:0); // 

    
    return x ;
}

Now we calculate for each Session where Re = {400,700,100} for each ER = {1.5,2.0,2.5} the friction coefficient and Nusselt number

var Reynolds = ((XNSEC_Control)myDb.Sessions[0].GetControl()).SelfDefinedHomotopyArray; // Obtain a list with Reynolds numbers used in the homotopy strategy
double[] ER = new double[]{2.0}; // All results are shown for a ER of two
double[] xCoordinates  = GenericBlas.Linspace(0,50, 1000); // points where values of solutions fields are going to be picked for calculation of the Nusselt numbers

// Results are stored for each (ER,Reynolds) solution pair
Dictionary<(double,double), double[]> cfs = new Dictionary<(double,double), double[]>();
Dictionary<(double,double), double[]> Nusselts = new Dictionary<(double,double), double[]>();
Dictionary<(double,double), double[]> ReattachmentPositions = new Dictionary<(double,double),  double[]>();
Dictionary<(double,double), double[]> ReattachmentPositionsTOP = new Dictionary<(double,double), double[]>();

foreach(double er in ER){
    var sess = myDb.Sessions.Where(s => s.Tags.Pick(0) == er.ToString()).SingleOrDefault(); // Session with ER = 2.0
    var control =  (XNSEC_Control)sess.GetControl();
    var ReynoldsArray = control.SelfDefinedHomotopyArray;
    Dictionary<int, double> ReynoldsDictionary = new Dictionary<int, double>();
    for(int i = 0; i <ReynoldsArray.Length;i++ ){
        ReynoldsDictionary.Add(i, ReynoldsArray[i]);
    }


    foreach(double Re in Reynolds){        
        int indexOfReynolds = ReynoldsDictionary.FirstOrDefault(x => x.Value == Re).Key;
        var timestep = sess.Timesteps[indexOfReynolds+1];

        // Calculation of the local Nusselt number          
        double tin = 10 +273; //Temperature inlet, K
        double thot = 40+273;//Temperature hot wall, K
        double tref = tin; // reference temperature
        double TIN = tin/tref; 
        double THOT = thot/tref; 
        var nuss = LocalNusseltNumber(timestep,xCoordinates, -(TIN-THOT), 1, 0 );
        Nusselts.Add((er,Re),nuss);

        //Calculation of the darcy friction factor
        double[] cf = LocalDarcyFactor(timestep, xCoordinates, 1,0,Re);
        cfs.Add((er,Re),cf);

        // // Reattachment positions 
        double[] X = GetPositionOfReattachment(timestep,1,0,er,top:false); 
        ReattachmentPositions.Add((er,Re),X);

        double[] XTop = GetPositionOfReattachment(timestep,1,0,er, top:true);
        ReattachmentPositionsTOP.Add((er,Re),XTop);

    }   
}

LineColors[] allColors = Enum.GetValues(typeof(LineColors)).Cast<LineColors>().ToArray();
PointTypes[] myPointTypes = new PointTypes[] { PointTypes.Diamond, PointTypes.Box, PointTypes.LowerTriangle, PointTypes.OpenLowerTriangle, };
LineColors[] myCollors = new LineColors[] { LineColors.Red, LineColors.Orange, LineColors.Blue };

First we compare our results with the ones given in the benchmark Biswas, G., M. Breuer, and F. Durst. “Backward-Facing Step Flows for Various Expansion Ratios at Low and Moderate Reynolds Numbers.” Journal of Fluids Engineering 126, no. 3 (July 12, 2004): 362–74. https://doi.org/10.1115/1.1760532.

In particular we compare results for $\text{Re} = 700$

Points of separation are obtained by finding the positions where the value of the gradient of velocityX normal to the wall changes its sign.

List<double> Rey = new List<double>();
List<double> Lenghths = new List<double>();

List<double> ReyTop = new List<double>();
List<double> Lenghths_left = new List<double>();
List<double> Lenghths_right = new List<double>();

foreach(var a in ReattachmentPositionsTOP){
    try{
        ReyTop.Add(a.Key.Item2);
        Lenghths_left.Add(a.Value[0]);
        Lenghths_right.Add(a.Value[1]);

        }catch(Exception e){
            Lenghths_left.Add(0);
            Lenghths_right.Add(0);
        continue;
        }
}
foreach(var a in ReattachmentPositions){
        double Reynolds = a.Key.Item2;
        double[] allLengths = a.Value;
        try{
            double primLength = allLengths[0];
            if(primLength < 0.3) //Try to filter values in order to obtain the length of reattachment after the aparition of the primary vortex  
                primLength = allLengths[1];                
            Rey.Add(Reynolds);
            Lenghths.Add(primLength);
        }catch(Exception e){
        continue;
        }
}

(14,26): warning CS0168: The variable 'e' is declared but never used



(29,26): warning CS0168: The variable 'e' is declared but never used

double nu = 15.52e-6;// Kinematic viscosity of air at 25°c, m2/s
double S = 15.0/1000; // Reference length according to the paper of XIA, meters
double vel = 700*nu/S; // mean velocity, m/s

double factor = 4*nu/(vel*vel/4.0);

Comparison of reattachment lengths for the primary vortex for ER =1.9423

var plot = new Plot2Ddata();
var fmt = new PlotFormat();
fmt.Style = Styles.Lines;

fmt.LineColor = myCollors[0];
plot.AddDataGroup("BoSSS", Rey.ToArray(), Lenghths.ToArray(), fmt);

// Reference data
MultidimensionalArray[] ReferenceDataNu = new MultidimensionalArray[1];

string path2 = Path.Combine(Directory.GetCurrentDirectory(), @"reatlengths_bot.txt"); // CF_RE400_ER1data for ER 1.94...
ReferenceDataNu[0] = IMatrixExtensions.LoadFromTextFile(path2);
var fmt4 = new PlotFormat();
fmt4.Style = Styles.Points;

fmt4.PointType = myPointTypes[0];// PointTypes.Diamond;
fmt4.PointSize = 0.5;
var xArr_nu = ReferenceDataNu[0].ExtractSubArrayShallow(-1, 0).To1DArray();
for(int i = 0; i < xArr_nu.Length; i++){
    xArr_nu[i]*=0.5; // REescale Reynolds for comparison with paper 
}

var yArr_nu = ReferenceDataNu[0].ExtractSubArrayShallow(-1, 1).To1DArray();
plot.AddDataGroup("reference",xArr_nu,yArr_nu,fmt4);
plot . PlotNow ()

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

var plot = new Plot2Ddata();
var fmt = new PlotFormat();
fmt.Style = Styles.Lines;

fmt.LineColor = myCollors[0];
plot.AddDataGroup("BoSSS_Deatach",ReyTop.ToArray(),Lenghths_left.ToArray(), fmt);
plot.AddDataGroup("BoSSS_ReAtach",ReyTop.ToArray(),Lenghths_right.ToArray(), fmt);

// Reference data
MultidimensionalArray[] ReferenceDataLeft = new MultidimensionalArray[1];
string path2 = Path.Combine(Directory.GetCurrentDirectory(), @"reatlengths_top_left.txt"); // CF_RE400_ER1data for ER 1.94...
ReferenceDataLeft[0] = IMatrixExtensions.LoadFromTextFile(path2);

MultidimensionalArray[] ReferenceDataRight = new MultidimensionalArray[1];
string path = Path.Combine(Directory.GetCurrentDirectory(), @"reatlengths_top_right.txt"); // CF_RE400_ER1data for ER 1.94...
ReferenceDataRight[0] = IMatrixExtensions.LoadFromTextFile(path);


var fmt4 = new PlotFormat();
fmt4.Style = Styles.Points;
fmt4.PointType = myPointTypes[0];// PointTypes.Diamond;
fmt4.PointSize = 0.5;


var xArr_left = ReferenceDataLeft[0].ExtractSubArrayShallow(-1, 0).To1DArray();
var yArr_left = ReferenceDataLeft[0].ExtractSubArrayShallow(-1, 1).To1DArray();
for(int i = 0; i < xArr_left.Length; i++){
    xArr_left[i]*=0.5; // REescale Reynolds for comparison with paper 
}


var xArr_right = ReferenceDataRight[0].ExtractSubArrayShallow(-1, 0).To1DArray();
var yArr_right = ReferenceDataRight[0].ExtractSubArrayShallow(-1, 1).To1DArray();
for(int i = 0; i < xArr_right.Length; i++){
    xArr_right[i]*=0.5; // REescale Reynolds for comparison with paper 
}



plot.AddDataGroup("reference",xArr_left,yArr_left,fmt4);
plot.AddDataGroup("reference",xArr_right,yArr_right,fmt4);
plot.PlotNow ()

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

double er = 2.0;
var sess = myDb.Sessions.Where(s => s.Tags.Pick(0) == er.ToString()).SingleOrDefault();

var control =  (XNSEC_Control)sess.GetControl();
var ReynoldsArray = control.SelfDefinedHomotopyArray;
Dictionary<int, double> ReynoldsDictionary = new Dictionary<int, double>();
for(int i = 0; i <ReynoldsArray.Length;i++ ){
    ReynoldsDictionary.Add(i, ReynoldsArray[i]);
}


double Re = 400;
int indexOfReynolds = ReynoldsDictionary.FirstOrDefault(x => x.Value == Re).Key;
var timestep = sess.Timesteps[indexOfReynolds+1];
var VelocityX = timestep.Fields.Where(f=> f.Identification == VariableNames.VelocityX).SingleOrDefault();
var VelocityY = timestep.Fields.Where(f=> f.Identification == VariableNames.VelocityY).SingleOrDefault();

double[] ys = GenericBlas.Linspace(0.0001,1.999,41);
double[] ys_real = GenericBlas.Linspace(0.0001,1.999,41);
double[] VelXArr7 = new double[ys.Length];
double[] VelYArr7 = new double[ys.Length];

double[] VelXArr15 = new double[ys.Length];
double[] VelYArr15 = new double[ys.Length];


for(int i = 0; i< ys.Length; i++){
    VelXArr7[i] = VelocityX.ProbeAt(new double[] {7.02*2,ys[i]})*0.5;
    VelXArr15[i] = VelocityX.ProbeAt(new double[] {15.0*2,ys[i]})*0.5;
    
    VelYArr7[i] = VelocityY.ProbeAt(new double[] {7*2,ys[i]})*0.5;
    VelYArr15[i] = VelocityY.ProbeAt(new double[] {15.0*2,ys[i]})*0.5;

    ys_real[i] = ys[i]/2.0 - 0.5;
}

var plot = new Plot2Ddata();
var fmt = new PlotFormat();
fmt.Style = Styles.Lines;

fmt.LineColor = myCollors[0];
plot.AddDataGroup("BoSSS",VelXArr7, ys_real,  fmt);
plot.AddDataGroup("BoSSS",VelXArr15, ys_real,  fmt);

// Reference data
MultidimensionalArray[] ReferenceData1 = new MultidimensionalArray[1];
MultidimensionalArray[] ReferenceData2 = new MultidimensionalArray[1];

string path1 = Path.Combine(Directory.GetCurrentDirectory(), @"Re400_u_x7_.txt"); // CF_RE400_ER1data for ER 1.94...
string path2 = Path.Combine(Directory.GetCurrentDirectory(), @"Re400_u_x15_.txt"); // CF_RE400_ER1data for ER 1.94...
ReferenceData1[0] = IMatrixExtensions.LoadFromTextFile(path1);
ReferenceData2[0] = IMatrixExtensions.LoadFromTextFile(path2);
var fmt4 = new PlotFormat();
fmt4.Style = Styles.Points;

fmt4.PointType = myPointTypes[0];// PointTypes.Diamond;
fmt4.PointSize = 0.5;
var xArr_1 = ReferenceData1[0].ExtractSubArrayShallow(-1, 0).To1DArray();
var xArr_2 = ReferenceData2[0].ExtractSubArrayShallow(-1, 0).To1DArray();
for(int i = 0; i < xArr_1.Length; i++){
    xArr_1[i]*=0.5; // REescale for comparison with paper  
}
for(int i = 0; i < xArr_2.Length; i++){
    xArr_2[i]*=0.5; // REescale for comparison with paper  
}

var yArr_1 = ReferenceData1[0].ExtractSubArrayShallow(-1, 1).To1DArray();
var yArr_2 = ReferenceData2[0].ExtractSubArrayShallow(-1, 1).To1DArray();
plot.AddDataGroup("reference",xArr_1,yArr_1,fmt4);
plot.AddDataGroup("reference",xArr_2,yArr_2,fmt4);
plot.PlotNow ()

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

var plot = new Plot2Ddata();
var fmt = new PlotFormat();
fmt.Style = Styles.Lines;

fmt.LineColor = myCollors[0];
// plot.AddDataGroup("BoSSS",VelYArr7, ys_real,  fmt);
plot.AddDataGroup("BoSSS",VelYArr15, ys_real,  fmt);

// Reference data
MultidimensionalArray[] ReferenceData1 = new MultidimensionalArray[1];
MultidimensionalArray[] ReferenceData2 = new MultidimensionalArray[1];

string path1 = Path.Combine(Directory.GetCurrentDirectory(), @"Re400_v_x7_.txt"); 
string path2 = Path.Combine(Directory.GetCurrentDirectory(), @"Re400_v_x15_.txt");
ReferenceData1[0] = IMatrixExtensions.LoadFromTextFile(path1);
ReferenceData2[0] = IMatrixExtensions.LoadFromTextFile(path2);
var fmt4 = new PlotFormat();
fmt4.Style = Styles.Points;

fmt4.PointType = myPointTypes[0];// PointTypes.Diamond;
fmt4.PointSize = 0.5;
var xArr_1 = ReferenceData1[0].ExtractSubArrayShallow(-1, 0).To1DArray();
var xArr_2 = ReferenceData2[0].ExtractSubArrayShallow(-1, 0).To1DArray();
for(int i = 0; i < xArr_1.Length; i++){
    xArr_1[i]*=0.5; // REescale for comparison with paper  
}
for(int i = 0; i < xArr_2.Length; i++){
    xArr_2[i]*=0.5; // REescale for comparison with paper  
}
var yArr_1 = ReferenceData1[0].ExtractSubArrayShallow(-1, 1).To1DArray();
var yArr_2 = ReferenceData2[0].ExtractSubArrayShallow(-1, 1).To1DArray();
// plot.AddDataGroup("reference",xArr_1,yArr_1,fmt4);
plot.AddDataGroup("reference",xArr_2,yArr_2,fmt4);
plot.PlotNow ()

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

string basepath = @"C:\tmp\NU_FD\";

private static void ToTxtFile(string fileName, double[] x, double[] y) {
    if(x.Length != y.Length)
        throw new Exception();
        using (var file = new StreamWriter(fileName)) {
        for(int i = 0; i< x.Length; i++){
            file.Write(x[i] + "\t" + y[i]);
            file.WriteLine();
        }
    }    
}

bool exportToTxt  = false;
if(exportToTxt){
    ToTxtFile(basepath + "BottomReatLength.txt", Rey.ToArray(), Lenghths.ToArray());
    ToTxtFile(basepath + "TopSepLength.txt",ReyTop.ToArray(),Lenghths_left.ToArray());
    ToTxtFile(basepath + "TopReatLength.txt", ReyTop.ToArray(),Lenghths_right.ToArray());

    ToTxtFile(basepath + "REF_BottomReatLength.txt",xArr_nu, yArr_nu);
    ToTxtFile(basepath + "REF_TopSepLength.txt",xArr_left, yArr_left);
    ToTxtFile(basepath + "REF_TopReatLength.txt",xArr_right, yArr_right);

    ToTxtFile(basepath + "VelXArr7.txt" ,VelXArr7, ys_real);
    ToTxtFile(basepath + "VelXArr15.txt" ,VelXArr15, ys_real);
    ToTxtFile(basepath + "VelXArr7_REF.txt" ,xArr_1, yArr_1);
    ToTxtFile(basepath + "VelXArr15_REF.txt" ,xArr_2, yArr_2);  
    

}

Plot2Ddata[,] PlotTableCfs = new Plot2Ddata[1, 3];
Plot2Ddata[,] PlotTableNusselts = new Plot2Ddata[1, 3];
int counter = 0;
double[] ReynoldsFromPaper = new double[]{700};
foreach (double Re in ReynoldsFromPaper) {
    var plotcf = new Plot2Ddata();
    var plotnuss = new Plot2Ddata();

    // BoSSS Data
    int cnt = 0;
    foreach (double er in ER) {
        var fmt = new PlotFormat();
        fmt.Style = Styles.Lines;

        fmt.LineColor = myCollors[cnt];
        plotcf.AddDataGroup("ER=" + er, xCoordinates, cfs[(er, Re)], fmt);
        plotnuss.AddDataGroup("ER=" + er, xCoordinates, Nusselts[(er, Re)], fmt);
       
        // Add reference data
        // Data avaliable for Re 400 700 1000 and ER 1.5 2.0 2.5
        // if(er == 1.5 || er == 2.0 || er == 2.5){
        MultidimensionalArray[] ReferenceDatacf = new MultidimensionalArray[1];
        MultidimensionalArray[] ReferenceDataNu = new MultidimensionalArray[1];
        string CurrentDocDir = Directory.GetCurrentDirectory();
        string erstr = ((int)(er * 10)).ToString();

        // string path1 = Path.Combine(CurrentDocDir, @"CF_RE" + Re + "_ER" + erstr + ".txt"); // CF_RE400_ER15 Data avaliable for Re 400 700 1000 and ER 1.5 2.0 2.5
        // string path2 = Path.Combine(CurrentDocDir, @"NU_RE" + Re + "_ER" + erstr + ".txt"); // CF_RE400_ER15
        string path1 = Path.Combine(CurrentDocDir, @"fd700.txt"); // CF_RE400_ER15 Data avaliable for Re 400 700 1000 and ER 1.5 2.0 2.5
        string path2 = Path.Combine(CurrentDocDir, @"NU__.txt"); // CF_RE400_ER15

        ReferenceDatacf[0] = IMatrixExtensions.LoadFromTextFile(path1);
        ReferenceDataNu[0] = IMatrixExtensions.LoadFromTextFile(path2);
        var fmt4 = new PlotFormat();
        fmt4.Style = Styles.Points;
        // fmt.Style = Styles.Lines;

        fmt4.PointType = myPointTypes[cnt];// PointTypes.Diamond;
        fmt4.PointSize = 0.5;
        

        var xArr_cf = ReferenceDatacf[0].ExtractSubArrayShallow(-1, 0).To1DArray();
        var yArr_cf = ReferenceDatacf[0].ExtractSubArrayShallow(-1, 1).To1DArray();
        var xArr_nu = ReferenceDataNu[0].ExtractSubArrayShallow(-1, 0).To1DArray();
        var yArr_nu = ReferenceDataNu[0].ExtractSubArrayShallow(-1, 1).To1DArray();
       


        plotcf.AddDataGroup(xArr_cf,yArr_cf,fmt4);
        plotcf.XrangeMin=xCoordinates.Min();
        plotcf.XrangeMax=xCoordinates.Max();

    
        plotnuss.AddDataGroup( xArr_nu,yArr_nu, fmt4);
        plotnuss.XrangeMin=xCoordinates.Min();
        plotnuss.XrangeMax=xCoordinates.Max();
    // }

        cnt++; 
    }
    PlotTableCfs[0, counter] = plotcf;
    PlotTableNusselts[0, counter] = plotnuss;
    counter++;
}

var gp = PlotTableCfs.ToGnuplot();
gp.PlotSVG(xRes:2000,yRes:400)

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

warning CS1701: Assuming assembly reference 'Microsoft.AspNetCore.Html.Abstractions, Version=2.2.0.0, Culture=neutral, PublicKeyToken=adb9793829ddae60' used by 'BoSSSpad' matches identity 'Microsoft.AspNetCore.Html.Abstractions, Version=6.0.0.0, Culture=neutral, PublicKeyToken=adb9793829ddae60' of 'Microsoft.AspNetCore.Html.Abstractions', you may need to supply runtime policy

var gp = PlotTableNusselts.ToGnuplot();
gp.PlotSVG(xRes:2000,yRes:400)

Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe

Note: In a Jupyter Worksheet, you must NOT have a trailing semicolon in order to see the plot on screen; otherwise, the output migth be surpressed.!

warning CS1701: Assuming assembly reference 'Microsoft.AspNetCore.Html.Abstractions, Version=2.2.0.0, Culture=neutral, PublicKeyToken=adb9793829ddae60' used by 'BoSSSpad' matches identity 'Microsoft.AspNetCore.Html.Abstractions, Version=6.0.0.0, Culture=neutral, PublicKeyToken=adb9793829ddae60' of 'Microsoft.AspNetCore.Html.Abstractions', you may need to supply runtime policy