Results published: Gutiérrez-Jorquera, Kummer: A fully coupled high-order discontinuous Galerkin method for diffusion flames in a low-Mach number framework
This example can be found in the source code repository as as HeatedCavity_RaSweepPostProc.ipynb.
One can directly load this into Jupyter to interactively work with the following code examples.
Note: First, BoSSS has to be loaded into the Jupyter kernel. Note:
In the following line, the reference to BoSSSpad.dll is required.
One must either set #r "BoSSSpad.dll" to something which is appropirate for the current computer
(e.g. C:\Program Files (x86)\FDY\BoSSS\bin\Release\net5.0\BoSSSpad.dll if working with the binary distribution),
or, if one is working with the source code, one must compile BoSSSpad
and put it side-by-side to this worksheet file
(from the original location in the repository, one can use the scripts getbossspad.sh, resp. getbossspad.bat).
#r "BoSSSpad.dll"
// #r "C:\BoSSS2\experimental\public\src\L4-application\BoSSSpad\bin\Release\net5.0\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();
BoSSSshell.WorkflowMgm.Init("HeatedCavity_RayleighSweepStudy");
Project name is set to 'HeatedCavity_RayleighSweepStudy'. Default Execution queue is chosen for the database. Opening existing database '\\fdygitrunner\ValidationTests\databases\HeatedCavity_RayleighSweepStudy'.
var sortedSessions = wmg.Sessions.OrderBy(sess => Convert.ToDouble(sess.KeysAndQueries["Rayleigh"]));
sortedSessions
#0: HeatedCavity_RayleighSweepStudy NaturalConvection_k5_DG3_Ra100 11/09/2023 13:17:51 082ec5d1... #1: HeatedCavity_RayleighSweepStudy NaturalConvection_k5_DG3_Ra1000 11/09/2023 13:18:08 cf5599f9... #2: HeatedCavity_RayleighSweepStudy NaturalConvection_k5_DG3_Ra10000 11/09/2023 13:18:27 b7690393... #3: HeatedCavity_RayleighSweepStudy NaturalConvection_k5_DG3_Ra100000 11/09/2023 13:18:47 543c6f35... #4: HeatedCavity_RayleighSweepStudy NaturalConvection_k5_DG3_Ra1000000 11/09/2023 13:19:06 9ceb8fb7... #5: HeatedCavity_RayleighSweepStudy NaturalConvection_k5_DG3_Ra10000000 11/09/2023 13:19:26 a6cc631b...
Dictionary<double,double> p0_solutions = new Dictionary<double,double>();
p0_solutions.Add(1e2, 0.9573);
p0_solutions.Add(1e3, 0.9381);
p0_solutions.Add(1e4, 0.9146);
p0_solutions.Add(1e5, 0.9220);
p0_solutions.Add(1e6, 0.9245);
p0_solutions.Add(1e7, 0.9226);
Dictionary<double,double> NusseltSolutions = new Dictionary<double,double>();
NusseltSolutions.Add(1e2, 0.9787);
NusseltSolutions.Add(1e3, 1.1077);
NusseltSolutions.Add(1e4, 2.2180);
NusseltSolutions.Add(1e5, 4.4800);
NusseltSolutions.Add(1e6, 8.6870);
NusseltSolutions.Add(1e7, 16.2400);
foreach(var sess in sortedSessions) {
double Rayleigh = Convert.ToDouble(sess.KeysAndQueries["Rayleigh"]);
Console.WriteLine("Rayleigh number:" + Rayleigh);
// Pick thermodynamic pressure value
var p0 = sess.Timesteps[1].Fields.Where(f => f.Identification == "ThermodynamicPressure").SingleOrDefault().ProbeAt(0.5,0.5);
double p0Reference = p0_solutions[Rayleigh];
double p0Error = Math.Abs(p0Reference - p0);
Console.WriteLine("Discrepancy of p0: "+ p0Error);
//NUnit.Framework.Assert.IsTrue(p0Error< 1e-2);
// Nusselt number comparison
double NusseltReference = NusseltSolutions[Rayleigh];
double leftvalue = Convert.ToDouble(sess.KeysAndQueries["NusseltNumber0"]) / 1.2;
double rightvalue = Convert.ToDouble(sess.KeysAndQueries["NusseltNumber1"]) / 1.2*(-1);
double residualValue = Convert.ToDouble(sess.KeysAndQueries["NusseltNumber2"]);
NUnit.Framework.Assert.IsTrue((Math.Abs(NusseltReference - leftvalue)) < 1.0); // big treshold for now...
NUnit.Framework.Assert.IsTrue((Math.Abs(NusseltReference - rightvalue)) < 1.0);
Console.WriteLine(Math.Abs(NusseltReference - leftvalue));
Console.WriteLine(Math.Abs(NusseltReference - rightvalue));
// Console.WriteLine("For Ra" + Rayleigh + " the error is" + Math.Abs(NusseltReference - rightvalue));
// // Console.WriteLine(Math.Abs(leftvalue));
// Console.WriteLine(Math.Abs(rightvalue));
}
Rayleigh number:100 Discrepancy of p0: 6.581895837509677E-05 2.4572760868535326E-05 7.965746323801426E-05 Rayleigh number:1000 Discrepancy of p0: 5.343703230820118E-05 4.064061141750841E-05 0.00013391039778398728 Rayleigh number:10000 Discrepancy of p0: 2.9823641442439097E-05 7.023894369284633E-05 0.002175898953845845 Rayleigh number:100000 Discrepancy of p0: 3.1333696267843436E-05 0.00040619267758756905 0.01906863956067273 Rayleigh number:1000000 Discrepancy of p0: 4.9440460990868296E-05 0.008351370214073839 0.12747154008564365 Rayleigh number:10000000 Discrepancy of p0: 0.00041519344518647916 0.09677531490992664 0.5822361199465593
public static class HelperMethods{
// Used for avoiding picking values outside the computational domain
public static double[] GetEpsilon(double xi, double yi, double xright, double ytop) {
double[] epsXY = new double[2];
double epsx = 0;
double epsy = 0;
double eps = 1e-7;
if (xi == xright && yi == ytop) {
epsx = -eps;
epsy = -eps;
} else if (yi == ytop) {
epsx = eps;
epsy = -eps;
} else if (xi == xright) {
epsx = -eps;
epsy = eps;
} else {
epsx = eps;
epsy = eps;
}
epsXY[0] = epsx;
epsXY[1] = epsy;
return epsXY;
}
}
double xleft = -0.5;
double xright = 0.5;
double ybot = -0.5;
double ytop = 0.5;
int nPoints = 60;
double eps = 1e-7;
double[] Rayleighs = new double[]{1e2,1e3,1e4,1e5,1e6,1e7};
double[] x_p = GenericBlas.Linspace(xleft,xright ,nPoints);
double[] y_p = GenericBlas.Linspace(ybot,ytop,3);
MultidimensionalArray[] ReferenceDataBot = new MultidimensionalArray[6];
MultidimensionalArray[] ReferenceDataMid = new MultidimensionalArray[6];
MultidimensionalArray[] ReferenceDataTop = new MultidimensionalArray[6];
string CurrentDocDir = Directory.GetCurrentDirectory();
for (int i = 0; i < 6; i++) {
string path1 = Path.Combine(CurrentDocDir, @"Temperatures_bot_Ra" + i + "_REF" + ".txt");
string path2 = Path.Combine(CurrentDocDir, @"Temperatures_mid_Ra" + i + "_REF" + ".txt");
string path3 = Path.Combine(CurrentDocDir, @"Temperatures_top_Ra" + i + "_REF" + ".txt");
ReferenceDataBot[i] = IMatrixExtensions.LoadFromTextFile(path1);
ReferenceDataMid[i] = IMatrixExtensions.LoadFromTextFile(path2);
ReferenceDataTop[i] = IMatrixExtensions.LoadFromTextFile(path3);
}
Plot2Ddata[,] PlotTable = new Plot2Ddata[2, 3];
int cnt2 = -1;
for (int iCol = 0; iCol < 2; iCol++) {
for (int iRow = 0; iRow < 3; iRow++) {
cnt2++;
double[] TemperatureBot = new double[x_p.Count()];
double[] TemperatureMid = new double[x_p.Count()];
double[] TemperatureTop = new double[x_p.Count()];
// var _fields = sortedSessions.Pick(cnt2).Timesteps.Pick(1).Fields;
var _fields = sortedSessions.Pick(cnt2).Timesteps.Pick(1).Fields;
var Temperature = _fields.Where(f => f.Identification == VariableNames.Temperature).SingleOrDefault();
double[] yS = new double[] { -0.5, 0.0, 0.5 };
for (int i = 0; i < x_p.Count(); i++) {
foreach (double yi in yS) {
double xi = x_p[i];
double[] epsXY = HelperMethods.GetEpsilon(xi, yi, xright, ytop);
var ti = Temperature.ProbeAt(xi + epsXY[0], yi + epsXY[1]);
if (yi == yS[0]) {
TemperatureBot[i] = ti;
} else if (yi == yS[1]) {
TemperatureMid[i] = ti;
} else if (yi == yS[2]) {
TemperatureTop[i] = ti;
} else {
throw new Exception();
}
}
}
// "Repair" x-coordinate to compare it with Benchmark
var xp2 = new double[x_p.Count()];
for(int i = 0; i< x_p.Count(); i++){
xp2[i] = x_p[i] +0.5;
}
LineColors[] allColors = Enum.GetValues(typeof(LineColors)).Cast<LineColors>().ToArray();
PointTypes[] myPointTypes = new PointTypes[]{ PointTypes.Diamond, PointTypes.Box, PointTypes.LowerTriangle,PointTypes.OpenLowerTriangle, };
var plot = new Plot2Ddata();
var fmt = new PlotFormat();
fmt.Style = Styles.Lines;
fmt.DashType = DashTypes.Solid;
fmt.LineColor = LineColors.Blue;
plot.AddDataGroup( "y=-0.5",xp2, TemperatureBot, fmt);
var fmt2 = new PlotFormat();
fmt2.Style = Styles.Lines;
fmt2.DashType = DashTypes.Solid;
fmt2.LineColor = LineColors.Green;
plot.AddDataGroup( "y=0.0",xp2, TemperatureMid, fmt2);
var fmt3 = new PlotFormat();
fmt3.Style = Styles.Lines;
fmt3.DashType = DashTypes.Solid;
fmt3.LineColor = LineColors.Red;
plot.AddDataGroup( "y=+0.5",xp2, TemperatureTop, fmt3);
var fmt4 = new PlotFormat();
fmt4.Style = Styles.Points;
fmt4.PointType = PointTypes.Diamond;
plot.AddDataGroup( "TemperatureBot",
ReferenceDataBot[cnt2].ExtractSubArrayShallow(-1, 0).To1DArray(),
ReferenceDataBot[cnt2].ExtractSubArrayShallow(-1, 1).To1DArray(),
fmt4);
var fmt5 = new PlotFormat();
fmt5.Style = Styles.Points;
fmt5.PointType = PointTypes.OpenDiamond;
plot.AddDataGroup( "TemperatureMid",
ReferenceDataMid[cnt2].ExtractSubArrayShallow(-1, 0).To1DArray(),
ReferenceDataMid[cnt2].ExtractSubArrayShallow(-1, 1).To1DArray(),
fmt5);
var fmt6 = new PlotFormat();
fmt6.Style = Styles.Points;
fmt6.PointType = PointTypes.UpperTriangle;
plot.AddDataGroup( "TemperatureTop",
ReferenceDataTop[cnt2].ExtractSubArrayShallow(-1, 0).To1DArray(),
ReferenceDataTop[cnt2].ExtractSubArrayShallow(-1, 1).To1DArray(),
fmt6);
plot.Title ="Ra:1e" + (cnt2 +2) ;
plot.TitleFont = 20;
// plot.XrangeMin = -0.5;
// plot.XrangeMax = 0.5;
plot.XrangeMin = 0.0;
plot.XrangeMax = 1.0;
plot.YrangeMin = 0.4;
plot.YrangeMax = 1.6;
plot.ShowLegend = cnt2 == 0;
plot.Xlabel = "x";
plot.Ylabel = "T";
// PlotTable[iRow, iCol] = plot;
PlotTable[iCol,iRow] = plot;
}
}
var gp = PlotTable.ToGnuplot();
gp.PlotSVG(xRes:1200,yRes:800)
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
double[] x_p2 = new double[]{0.0};
double[] y_p2 = GenericBlas.Linspace(ybot,ytop,nPoints);
MultidimensionalArray[] ReferenceDataVelocityX = new MultidimensionalArray[6];
string CurrentDocDir = Directory.GetCurrentDirectory();
for (int i = 0; i < 6; i++) {
string path1 = String.Concat(CurrentDocDir, @"\VelocityX" + i + "_REF", ".txt");
ReferenceDataVelocityX[i] = IMatrixExtensions.LoadFromTextFile(path1);
}
Plot2Ddata[,] PlotTable = new Plot2Ddata[2, 3];
int cnt2 = -1;
for (int iCol = 0; iCol < 2; iCol++) {
for (int iRow = 0; iRow < 3; iRow++) {
cnt2++;
double[] VelocityX = new double[y_p2.Count()];
var _fields = sortedSessions.Pick(cnt2).Timesteps.Pick(1).Fields;
var VelocityXField = _fields.Where(f => f.Identification == VariableNames.VelocityX).SingleOrDefault();
for (int i = 0; i < y_p2.Count(); i++) {
double xi = 0.0; // centerline
double yi = y_p2[i];
double[] epsXY = HelperMethods.GetEpsilon(xi, yi, xright, ytop);
var ti = VelocityXField.ProbeAt(xi + epsXY[0], yi + epsXY[1]);
VelocityX[i] = ti;
}
// "Repair" x-coordinate to compare it with Benchmark
var yp3 = new double[y_p2.Count()];
for(int i = 0; i< y_p2.Count(); i++){
yp3[i] = y_p2[i] +0.5;
}
LineColors[] allColors = Enum.GetValues(typeof(LineColors)).Cast<LineColors>().ToArray();
PointTypes[] myPointTypes = new PointTypes[]{ PointTypes.Diamond, PointTypes.Box, PointTypes.LowerTriangle,PointTypes.OpenLowerTriangle, };
var plot = new Plot2Ddata();
var fmt = new PlotFormat();
fmt.Style = Styles.Lines;
fmt.DashType = DashTypes.Solid;
fmt.LineColor = LineColors.Blue;
plot.AddDataGroup( "x=0.0",yp3, VelocityX, fmt);
var fmt4 = new PlotFormat();
fmt4.Style = Styles.Points;
fmt4.PointType = PointTypes.Diamond;
plot.AddDataGroup( "VelocityX",
ReferenceDataVelocityX[cnt2].ExtractSubArrayShallow(-1, 0).To1DArray(),
ReferenceDataVelocityX[cnt2].ExtractSubArrayShallow(-1, 1).To1DArray(),
fmt4);
plot.Title ="Ra:1e" + (cnt2 +2) ;
plot.TitleFont = 20;
// plot.XrangeMin = -0.5;
// plot.XrangeMax = 0.5;
plot.XrangeMin = 0.0;
plot.XrangeMax = 1.0;
plot.YrangeMin = -0.4;
plot.YrangeMax = 0.4;
plot.ShowLegend = cnt2 == 0;
plot.Xlabel = "y";
plot.Ylabel = "u_x";
// PlotTable[iRow, iCol] = plot;
PlotTable[iCol,iRow] = plot;
}
}
var gp = PlotTable.ToGnuplot();
gp.PlotSVG(xRes:1200,yRes:800)
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
double[] x_p2 = GenericBlas.Linspace(xleft,xright,nPoints);
double[] y_p2 = new double[]{0.0};
MultidimensionalArray[] ReferenceDataVelocityY = new MultidimensionalArray[6];
string CurrentDocDir = Directory.GetCurrentDirectory();
for (int i = 0; i < 6; i++) {
string path1 = String.Concat(CurrentDocDir, @"\VelocityY" + i + "_REF", ".txt");
ReferenceDataVelocityY[i] = IMatrixExtensions.LoadFromTextFile(path1);
}
Plot2Ddata[,] PlotTable = new Plot2Ddata[2, 3];
int cnt2 = -1;
for (int iCol = 0; iCol < 2; iCol++) {
for (int iRow = 0; iRow < 3; iRow++) {
cnt2++;
double[] VelocityY = new double[x_p2.Count()];
var _fields = sortedSessions.Pick(cnt2).Timesteps.Pick(1).Fields;
var VelocityYField = _fields.Where(f => f.Identification == VariableNames.VelocityY).SingleOrDefault();
for (int i = 0; i < x_p2.Count(); i++) {
double xi = x_p2[i]; // centerline
double yi = 0.0; // centerline
double[] epsXY = HelperMethods.GetEpsilon(xi, yi, xright, ytop);
var ti = VelocityYField.ProbeAt(xi + epsXY[0], yi + epsXY[1]);
VelocityY[i] = ti;
}
// "Repair" x-coordinate to compare it with Benchmark
var xp3 = new double[x_p2.Count()];
for(int i = 0; i< x_p2.Count(); i++){
xp3[i] = x_p2[i] +0.5;
}
LineColors[] allColors = Enum.GetValues(typeof(LineColors)).Cast<LineColors>().ToArray();
PointTypes[] myPointTypes = new PointTypes[]{ PointTypes.Diamond, PointTypes.Box, PointTypes.LowerTriangle,PointTypes.OpenLowerTriangle, };
var plot = new Plot2Ddata();
var fmt = new PlotFormat();
fmt.Style = Styles.Lines;
fmt.DashType = DashTypes.Solid;
fmt.LineColor = LineColors.Blue;
plot.AddDataGroup( "y=0.0",xp3, VelocityY, fmt);
var fmt4 = new PlotFormat();
fmt4.Style = Styles.Points;
fmt4.PointType = PointTypes.Diamond;
plot.AddDataGroup( "VelocityY",
ReferenceDataVelocityY[cnt2].ExtractSubArrayShallow(-1, 0).To1DArray(),
ReferenceDataVelocityY[cnt2].ExtractSubArrayShallow(-1, 1).To1DArray(),
fmt4);
plot.Title ="Ra:1e" + (cnt2 +2) ;
plot.TitleFont = 20;
// plot.XrangeMin = -0.5;
// plot.XrangeMax = 0.5;
plot.XrangeMin = 0.0;
plot.XrangeMax = 1.0;
plot.YrangeMin = -0.4;
plot.YrangeMax = 0.4;
plot.ShowLegend = cnt2 == 0;
plot.Xlabel = "x";
plot.Ylabel = "u_y";
// PlotTable[iRow, iCol] = plot;
PlotTable[iCol,iRow] = plot;
}
}
var gp = PlotTable.ToGnuplot();
gp.PlotSVG(xRes:1200,yRes:800)
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