No BoSSS specific prerequisites are needed to complete this tutorial.
First, we define two functions: $g_1$ is continuous, $g_2$ has a discontinuity at $x = \pi$ in the first derivative
$$ g_1(x) := \sin(x), $$$$g_2(x) := \vert \sin(x) \vert ,$$The function argument is a vector $x \in R^n$, consisting only of one entry since we are working in a one dimensional space.
BoSSS however supports 1D, 2D and 3D, so the spatial coordinate is a general vector.
Note:
ue2Basics.ipynb
.
One can directly load this into Jupyter to interactively work with the following code examples.BoSSSpad.dll
is required.
You must either set #r "BoSSSpad.dll"
to something which is appropirate for your computer
(e.g. C:\Program Files (x86)\FDY\BoSSS\bin\Release\net5.0\BoSSSpad.dll
if you installed the binary distribution),
or, if you are working with the source code, you must compile BoSSSpad
and put it side-by-side to this worksheet file
(from the original location in the repository, you can use the scripts getbossspad.sh
, resp. getbossspad.bat
).#r "BoSSSpad.dll"
using System;
using System.Collections.Generic;
using System.Linq;
using ilPSP;
using ilPSP.Utils;
using BoSSS.Platform;
using BoSSS.Platform.LinAlg;
using BoSSS.Foundation;
using BoSSS.Foundation.XDG;
using BoSSS.Foundation.Grid;
using BoSSS.Foundation.Grid.Classic;
using BoSSS.Foundation.Grid.RefElements;
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.AdvancedSolvers;
using BoSSS.Solution.Gnuplot;
using BoSSS.Application.BoSSSpad;
using BoSSS.Application.XNSE_Solver;
using static BoSSS.Application.BoSSSpad.BoSSSshell;
BoSSS.Application.BoSSSpad.BoSSSshell.Init();
using NUnit.Framework;
First, we plot the functions that are defined above over the interval $(0 ,2 \pi)$ with 1000 sampling points using a Gnuplot-object.
Func<double[],double> g1 = (X => Math.Sin(X[0]));
// continuous, smooth
Func<double[],double> g2 = (X => Math.Abs(Math.Sin(X[0])));
// continuous, non-smooth
We define equidistant sampling points...
double[] x = GenericBlas.Linspace(0, 2.0*Math.PI, 1000);
...and compute the function values. In the loop, we have to convert the scalar x[i] into an array with one element, since $g_1$ has to be feed with arrays.
double[] g1_values = new double[x.Length];
for(int i = 0; i < x.Length; i++) {
g1_values[i] = g1(new[] { x[i]} );
}
/// Instead of loops, we can also use Linq-functions:
double[] g2_values = x.Select(x => g2(new []{ x })).ToArray();
For now, we are using the simple plotting interface, which supports Matlab-Style format specifiers and color names.
(More advanced plots can be produced with Plot2Ddata and/or Gnuplot classes)
Note¶
In order to obtain an output for the plot or any other command, there must not be a semicolon ; at the end of the line!
Plot(X1:x, Y1:g1_values, Name1:"function g1", Format1:"--red",
X2:x, Y2:g2_values, Name2:"function g2", Format2:"-.blue")
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.!
Next, we create a grid which has a cell boundary exactly at the position of the discontinuity of $g_2$.
var Nodes1 = new double[] {0, 2, Math.PI, 4.5, 2*Math.PI };
var Grid1 = Grid1D.LineGrid(Nodes1);
We can get the total number of cells by using the following command:
Grid1.NumberOfCells
The recently created grid-object is not directly usable because it contains only the nodes of the grid.
We have to create a GridData-object which provides all necessary transformation metrics, etc. .
var gdata1 = new GridData(Grid1);
At this point, we are able to create the so-called DG fields to approximate $g_1$ on grid1.
Therefore, we project $g_1$ onto grid1 using polynomial orders of $p=2$ and $p=8$.
var g1_grid1_p2 = new SinglePhaseField(new Basis(gdata1, 2), "g1 with p2 at Grid 1");
g1_grid1_p2.ProjectField(g1);
var g1_grid1_p8 = new SinglePhaseField(new Basis(gdata1, 8), "g1 with p8 at Grid 1");
g1_grid1_p8.ProjectField(g1);
Now, let us plot the projected solution for $p=2$.
By using the upsampling parameter, we can determine the amount of sampling points per cell.
var upsampling = 20;
var gp1 = new Gnuplot();
Using gnuplot: C:\Program Files (x86)\FDY\BoSSS\bin\native\win\gnuplot-gp510-20160418-win32-mingw\gnuplot\bin\gnuplot.exe
gp1.PlotField(g1_grid1_p2,
new PlotFormat(lineColor: (LineColors)(1)),
upsampling);
gp1.PlotNow() // shows the plot
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.!
Next, we learn how to compute the $L^2$-error for both approximations of $g_1$ with different polynomial degrees:
g1_grid1_p2.L2Error(g1)
g1_grid1_p8.L2Error(g1)
Now, we plot the point-wise error for the approximation of $g_1$ on grid1 with a polynomial degree of 8.
int K = 20; // number of points per cell
var gp2 = new Gnuplot();
gp2.PlotLogError(g1_grid1_p8, g1, "g1 with p8 at Grid 1", 20,
new PlotFormat(lineColor: (LineColors)(1)));
gp2.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.!
We investigate the decay behavior of the DG modes for smooth and non-smooth functions.
For this purpose, we create a second grid which has the discontinuity of $g_2$ within a cell and project $g_2$ onto this grid like mentioned above.
var Nodes2 = new double[] {0, 2, 4.5, 2*Math.PI };
var Grid2 = Grid1D.LineGrid(Nodes2);
var gdata2 = new GridData(Grid2);
var g2_grid2_p8 = new SinglePhaseField(new Basis(gdata2, 8), "g2_p8 at Grid2");
g2_grid2_p8.ProjectField(g2);
The cell coordinates can be extracted by using the Coordinates parameter.
double[] cell1 = g2_grid2_p8.Coordinates.GetRow(1);
// coord. in cell 1 (with kink)
double[] cell0 = g2_grid2_p8.Coordinates.GetRow(0);
// coord. in cell 0 (smooth)
double[] cell2 = g2_grid2_p8.Coordinates.GetRow(2);
// coord. in cell 2 (smooth)
Only the absolute value shall be plotted. We use a for-loop to replace the data in cell0, cell1 and cell2 by their absolute values.
for(int i = 0; i < cell0.Length; i++) {
cell0[i] = Math.Abs(cell0[i]);
cell1[i] = Math.Abs(cell1[i]);
cell2[i] = Math.Abs(cell2[i]);
}
Plot(X1:null, Y1:cell1, Name1:"disc. cell", Format1:"*-magenta",
X2:null, Y2:cell0, Name2:"cell0", Format2:"o-blue",
X3:null, Y3:cell2, Name3:"cell2", Format3:"o-red")
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.!
Using a shortcut for the for-loop above, the absolute values in cell0 can also be stored using the following command:
double[] cell0 = g2_grid2_p8.Coordinates.GetRow(0)
.Select(d => Math.Abs(d)).ToArray();
Now, we would like to plot the logarithm (use Math.Log10(...)) of the absolute values of the DG coordinates.
Plot(X1:null, Y1:cell1, Name1:"disc. cell", Format1:"*-magenta",
X2:null, Y2:cell0, Name2:"cell0", Format2:"o-blue",
X3:null, Y3:cell2, Name3:"cell2", Format3:"o-red",
logY:true)
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.!
In this section, we learn how to perform a convergence study for $g_2$ for two different sequences of grid resolutions and different polynomial orders.
Therefore, we define two different sequences of grid resolutions:
int[][] ResSeq = new int[2][];
Grid resolutions so that the kink in g2 is located at a cell boundary:
ResSeq[0] = new int[] { 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 };
Grid resolutions so that the kink in g2 is located within a cell:
ResSeq[1] = new int[] { 3, 7, 15, 31, 63, 127, 255, 511, 1023, 2047 };
We save our errors into a multidimensional array by looping over
var Errors = MultidimensionalArray.Create(2, 5, ResSeq[0].Length);
for(int i = 0; i < 2; i++) { // loop over the resolution sequence
for(int p = 0; p <= 4; p++) { // loop over polynomial orders
for(int k = 0; k < ResSeq[i].Length; k++) { // loop over different resolutions
Console.Write("polynomial order {1}"+
",\tResolution {0}... ", ResSeq[i][k], p);
var grid = Grid1D.LineGrid(GenericBlas.Linspace(0,
2.0*Math.PI, ResSeq[i][k] + 1));
// number of nodes == number of cells + 1
var gData = new GridData(grid);
var g2_h = new SinglePhaseField(new Basis(gData, p));
g2_h.ProjectField(g2);
Errors[i,p,k] = g2_h.L2Error(g2);
Console.WriteLine("\tdone: L2 error is {0:0.###e-00}.", Errors[i,p,k]);
}
}
}
polynomial order 0, Resolution 4... done: L2 error is 1.791e-01. polynomial order 0, Resolution 8... done: L2 error is 4.536e-02. polynomial order 0, Resolution 16... done: L2 error is 1.138e-02. polynomial order 0, Resolution 32... done: L2 error is 2.846e-03. polynomial order 0, Resolution 64... done: L2 error is 7.118e-04. polynomial order 0, Resolution 128... done: L2 error is 1.779e-04. polynomial order 0, Resolution 256... done: L2 error is 4.449e-05. polynomial order 0, Resolution 512... done: L2 error is 1.112e-05. polynomial order 0, Resolution 1024... done: L2 error is 2.781e-06. polynomial order 0, Resolution 2048... done: L2 error is 6.951e-07. polynomial order 1, Resolution 4... done: L2 error is 2.155e-02. polynomial order 1, Resolution 8... done: L2 error is 2.739e-03. polynomial order 1, Resolution 16... done: L2 error is 3.438e-04. polynomial order 1, Resolution 32... done: L2 error is 4.302e-05. polynomial order 1, Resolution 64... done: L2 error is 5.379e-06. polynomial order 1, Resolution 128... done: L2 error is 6.724e-07. polynomial order 1, Resolution 256... done: L2 error is 8.405e-08. polynomial order 1, Resolution 512... done: L2 error is 1.051e-08. polynomial order 1, Resolution 1024... done: L2 error is 1.313e-09. polynomial order 1, Resolution 2048... done: L2 error is 1.642e-10. polynomial order 2, Resolution 4... done: L2 error is 2.112e-03. polynomial order 2, Resolution 8... done: L2 error is 1.34e-04. polynomial order 2, Resolution 16... done: L2 error is 8.405e-06. polynomial order 2, Resolution 32... done: L2 error is 5.258e-07. polynomial order 2, Resolution 64... done: L2 error is 3.287e-08. polynomial order 2, Resolution 128... done: L2 error is 2.055e-09. polynomial order 2, Resolution 256... done: L2 error is 1.284e-10. polynomial order 2, Resolution 512... done: L2 error is 8.027e-12. polynomial order 2, Resolution 1024... done: L2 error is 5.029e-13. polynomial order 2, Resolution 2048... done: L2 error is 3.274e-14. polynomial order 3, Resolution 4... done: L2 error is 1.665e-04. polynomial order 3, Resolution 8... done: L2 error is 5.273e-06. polynomial order 3, Resolution 16... done: L2 error is 1.653e-07. polynomial order 3, Resolution 32... done: L2 error is 5.171e-09. polynomial order 3, Resolution 64... done: L2 error is 1.616e-10. polynomial order 3, Resolution 128... done: L2 error is 5.051e-12. polynomial order 3, Resolution 256... done: L2 error is 1.576e-13. polynomial order 3, Resolution 512... done: L2 error is 4.887e-15. polynomial order 3, Resolution 1024... done: L2 error is 1.217e-15. polynomial order 3, Resolution 2048... done: L2 error is 1.008e-15. polynomial order 4, Resolution 4... done: L2 error is 1.094e-05. polynomial order 4, Resolution 8... done: L2 error is 1.73e-07. polynomial order 4, Resolution 16... done: L2 error is 2.711e-09. polynomial order 4, Resolution 32... done: L2 error is 4.24e-11. polynomial order 4, Resolution 64... done: L2 error is 6.647e-13. polynomial order 4, Resolution 128... done: L2 error is 1.267e-14. polynomial order 4, Resolution 256... done: L2 error is 3.253e-15. polynomial order 4, Resolution 512... done: L2 error is 3.138e-15. polynomial order 4, Resolution 1024... done: L2 error is 3.193e-15. polynomial order 4, Resolution 2048... done: L2 error is 3.098e-15. polynomial order 0, Resolution 3... done: L2 error is 8.806e-01. polynomial order 0, Resolution 7... done: L2 error is 2.499e-01. polynomial order 0, Resolution 15... done: L2 error is 7.914e-02. polynomial order 0, Resolution 31... done: L2 error is 2.65e-02. polynomial order 0, Resolution 63... done: L2 error is 9.121e-03. polynomial order 0, Resolution 127... done: L2 error is 3.182e-03. polynomial order 0, Resolution 255... done: L2 error is 1.117e-03. polynomial order 0, Resolution 511... done: L2 error is 3.938e-04. polynomial order 0, Resolution 1023... done: L2 error is 1.39e-04. polynomial order 0, Resolution 2047... done: L2 error is 4.91e-05. polynomial order 1, Resolution 3... done: L2 error is 2.415e-01. polynomial order 1, Resolution 7... done: L2 error is 6.353e-02. polynomial order 1, Resolution 15... done: L2 error is 2e-02. polynomial order 1, Resolution 31... done: L2 error is 6.713e-03. polynomial order 1, Resolution 63... done: L2 error is 2.316e-03. polynomial order 1, Resolution 127... done: L2 error is 8.09e-04. polynomial order 1, Resolution 255... done: L2 error is 2.843e-04. polynomial order 1, Resolution 511... done: L2 error is 1.002e-04. polynomial order 1, Resolution 1023... done: L2 error is 3.538e-05. polynomial order 1, Resolution 2047... done: L2 error is 1.25e-05. polynomial order 2, Resolution 3... done: L2 error is 2.622e-01. polynomial order 2, Resolution 7... done: L2 error is 7.053e-02. polynomial order 2, Resolution 15... done: L2 error is 2.23e-02. polynomial order 2, Resolution 31... done: L2 error is 7.494e-03. polynomial order 2, Resolution 63... done: L2 error is 2.586e-03. polynomial order 2, Resolution 127... done: L2 error is 9.033e-04. polynomial order 2, Resolution 255... done: L2 error is 3.175e-04. polynomial order 2, Resolution 511... done: L2 error is 1.119e-04. polynomial order 2, Resolution 1023... done: L2 error is 3.951e-05. polynomial order 2, Resolution 2047... done: L2 error is 1.396e-05. polynomial order 3, Resolution 3... done: L2 error is 1.263e-01. polynomial order 3, Resolution 7... done: L2 error is 3.488e-02. polynomial order 3, Resolution 15... done: L2 error is 1.109e-02. polynomial order 3, Resolution 31... done: L2 error is 3.731e-03. polynomial order 3, Resolution 63... done: L2 error is 1.288e-03. polynomial order 3, Resolution 127... done: L2 error is 4.499e-04. polynomial order 3, Resolution 255... done: L2 error is 1.581e-04. polynomial order 3, Resolution 511... done: L2 error is 5.574e-05. polynomial order 3, Resolution 1023... done: L2 error is 1.968e-05. polynomial order 3, Resolution 2047... done: L2 error is 6.952e-06. polynomial order 4, Resolution 3... done: L2 error is 1.338e-01. polynomial order 4, Resolution 7... done: L2 error is 3.703e-02. polynomial order 4, Resolution 15... done: L2 error is 1.178e-02. polynomial order 4, Resolution 31... done: L2 error is 3.963e-03. polynomial order 4, Resolution 63... done: L2 error is 1.368e-03. polynomial order 4, Resolution 127... done: L2 error is 4.778e-04. polynomial order 4, Resolution 255... done: L2 error is 1.679e-04. polynomial order 4, Resolution 511... done: L2 error is 5.92e-05. polynomial order 4, Resolution 1023... done: L2 error is 2.09e-05. polynomial order 4, Resolution 2047... done: L2 error is 7.384e-06.
/// NUnit test (few random tests)
Assert.LessOrEqual(Errors[1,4,9],8E-06);
Assert.LessOrEqual(Errors[1,3,9],7.5E-06);
Assert.LessOrEqual(Errors[1,2,9],2E-05);
Assert.LessOrEqual(Errors[1,1,9],2E-05);
Assert.LessOrEqual(Errors[0,3,9],1E-12);
Assert.LessOrEqual(Errors[0,3,9],1E-12);
Assert.LessOrEqual(Errors[1,4,0],0.25);
Assert.LessOrEqual(Errors[0,0,0],0.2);
Assert.LessOrEqual(Errors[0,3,0],1E-03);
We plot the error for the grids which have the kink at the cell boundary, there we reach spectral convergence:
var hValues = ResSeq[0].Select(J => Math.PI*2.0/J);
Plot(X1:hValues, Y1:Errors.ExtractSubArrayShallow(0,0,-1).To1DArray(),
Name1:"grid1,p0", Format1:"o-red",
X2:hValues, Y2:Errors.ExtractSubArrayShallow(0,1,-1).To1DArray(),
Name2:"grid1,p1", Format2:"o-blue",
X3:hValues, Y3:Errors.ExtractSubArrayShallow(0,2,-1).To1DArray(),
Name3:"grid1,p2", Format3:"o-green",
X4:hValues, Y4:Errors.ExtractSubArrayShallow(0,3,-1).To1DArray(),
Name4:"grid1,p3", Format4:"o-magenta",
X5:hValues, Y5:Errors.ExtractSubArrayShallow(0,4,-1).To1DArray(),
Name5:"grid1,p4", Format5:"o-orange",
logX:true, logY:true)
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.!
Now we plot the error for the grids which have the kink within a cell; due to the low regularity, the convergence of the DG method degenerates:
var hValues = ResSeq[0].Select(J => Math.PI*2.0/J);
Plot(X1:hValues, Y1:Errors.ExtractSubArrayShallow(1,0,-1).To1DArray(),
Name1:"grid1,p0", Format1:"o-red",
X2:hValues, Y2:Errors.ExtractSubArrayShallow(1,1,-1).To1DArray(),
Name2:"grid1,p1", Format2:"o-blue",
X3:hValues, Y3:Errors.ExtractSubArrayShallow(1,2,-1).To1DArray(),
Name3:"grid1,p2", Format3:"o-green",
X4:hValues, Y4:Errors.ExtractSubArrayShallow(1,3,-1).To1DArray(),
Name4:"grid1,p3", Format4:"o-magenta",
X5:hValues, Y5:Errors.ExtractSubArrayShallow(1,4,-1).To1DArray(),
Name5:"grid1,p4", Format5:"o-orange",
logX:true, logY:true)
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.!
This tutorial addressed the very basics of setting up a BoSSS~application, namely grid instantiation, the $L^2$-projection of functions onto the DG space and performing a spatial convergence study. Where do you go from here? We recommend that you continue with other relevant basics as provided in the tutorials dealing with the creation of a spatial operator, explicit time integration and the implementation of numerical fluxes.