Introduction

#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;
Init();

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

Note:

• Setting Boundary values and initial values is similar;
• For most solvers, inital and boundary values are set the same way;
• We will use the incompressible solver as an example:

using BoSSS.Application.XNSE_Solver;

Create a control object:

var C = new XNSE_Control();

1 From Formulas

If the Formula is simple enough to be represented by C# code, it can be embedded in the control file.

However, the code bust be put into a string, since it is not possible to serialize classes/objects from the notebook into a control object:

string code =  
 "static class MyInitialValue {" // class must be static!
 
    // Warning: static constants are allowed,
    // but any changes outside of the current text box in BoSSSpad
    // will not be recorded for the code that is passed to the solver.
+ " public static double alpha = 0.7;"
 
    // a method, which should be used for an initial value,
    // must be static!
+ " public static double VelocityX(double[] X, double t) {"
+ "    double x = X[0];"
+ "    double y = X[1];"
+ "   return Math.Sin(x*y*alpha);"
+ " }"  
+ "}";

var fo = new BoSSS.Solution.Control.Formula("MyInitialValue.VelocityX", 
                                             true, code);

Use the BoSSSpad-intrinsic GetFormulaObject to set tie inital value:

C.AddInitialValue("VelocityX", fo);

/// Deprecated:
/// Note: such a declaration is very restrictive;
/// \code{GetFormulaObject} works only for 
/// \begin{itemize}
/// \item a static class
/// \item no dependence on any external parameters
/// \end{itemize}
/// E.g. the following code would only change the behavior in BoSSSpad,
/// but not the code that is passed to the solver:

//Deprecated:
//MyInitialValue.alpha = 0.5;
//MyInitialValue.VelocityX(new double[]{ 0.5, 0.5 }, 0.0);

C.InitialValues["VelocityX"].Evaluate(new double[]{ 0.5, 0.5 }, 0.0)

0.17410813759359595

2 Advanced functions

Some more advanced mathematical functions, e.g. Jacobian elliptic functions $\text{sn}(u|m)$, $\text{cn}(u|m)$ and $\text{dn}(u|m)$ are available throug the GNU Scientific Library, for which BoSSS provides bindings, see e.g. BoSSS.Platform.GSL.gsl_sf_elljac_e

2.1 From MATLAB code

Asssume e.g. the following MATLAB code; obviously, this could
also be implemented in C#, we yust use something smple for demonstration:

string[] MatlabCode = new string[] {
@"[n,d2] =  size(X_values);",
@"u=zeros(2,n);",
@"for k=1:n",
@"X=[X_values(k,1),X_values(k,2)];",
@"",
@"u_x_main     = -(-sqrt(X(1).^ 2 + X(2).^ 2) / 0.3e1 + 0.4e1 / 0.3e1 * (X(1).^ 2 + X(2).^ 2) ^ (-0.1e1 / 0.2e1)) * sin(atan2(X(2), X(1)));",
@"u_y_main     = (-sqrt(X(1).^ 2 + X(2).^ 2) / 0.3e1 + 0.4e1 / 0.3e1 * (X(1).^ 2 + X(2).^ 2) ^ (-0.1e1 / 0.2e1)) * cos(atan2(X(2), X(1)));",
@"",   
@"u(1,k)=u_x_main;",
@"u(2,k)=u_y_main;",
@"end" };

We can evaluate this code in BoSSS using the MATLAB connector; We encapsulate it in a ScalarFunction which allows vectorized evaluation (multiple evaluatiuons in one function call) e of some function.

This is much more efficient, since there will be significant overhead for calling MATLAB (starting MATLAB, checking the license, transfering data, etc.).

using ilPSP.Connectors.Matlab;

ScalarFunction VelocityXInitial = 
delegate(MultidimensionalArray input, MultidimensionalArray output) {
    int N          = input.GetLength(0); // number of points which we evaluate 
    //                                      at once.
    var output_vec = MultidimensionalArray.Create(2, N); // the MATLAB code
    //                        returns an entire vector.
    using(var bmc = new BatchmodeConnector()) {
       bmc.PutMatrix(input,"X_values");
 
       foreach(var line in MatlabCode) {
           bmc.Cmd(line);   
       }
 
       bmc.GetMatrix(output_vec, "u");
 
       bmc.Execute(); // Note: 'Execute' has to be *after* 'GetMatrix'
   }
   output.Set(output_vec.ExtractSubArrayShallow(0,-1)); // extract row 0 from 
   //                       'output_vec' and store it in 'output'
};

We test our implementation:

var inputTest = MultidimensionalArray.Create(3,2); // set some test values for input
inputTest.SetColumn(0, GenericBlas.Linspace(1,2,3));
inputTest.SetColumn(1, GenericBlas.Linspace(2,3,3));
 
var outputTest = MultidimensionalArray.Create(3); // allocate memory for output

VelocityXInitial(inputTest, outputTest);

We recive the following velocity values for our input coordinates:

outputTest.To1DArray()

index	value
0	0.1333333333333334
1	0.4411764705882353
2	0.6923076923076924

Projecting the MATLAB function to a DG field

As for a standard calculation, we create a mesh, save it to some database and set the mesh in the control object.

var nodes        = GenericBlas.Linspace(1,2,11);
GridCommons grid = Grid2D.Cartesian2DGrid(nodes,nodes);

var db = CreateTempDatabase();

Creating database 'C:\Users\jenkinsci\AppData\Local\Temp\1398882499'.

db.SaveGrid(ref grid);

C.SetGrid(grid);

We create a DG field for the $x$-velocity on our grid:

var gdata = new GridData(grid);

var b = new Basis(gdata, 3); // use DG degree 2

var VelX = new SinglePhaseField(b,"VelocityX"); // important: name the DG field
//                                 equal to initial value name

Finally, we are able to project the MATLAB function onto the DG field:

//VelX.ProjectField(VelocityXInitial);

One might want to check the data visually, so it can be exported in the usual fashion

//Tecplot("initial",0.0,2,VelX);

Storing the initial value in the database and linking it in the control object

The DG field with the initial value can be stored in the database. this will create a dummy session.

BoSSSshell.WorkflowMgm.Init("TestProject");

Project name is set to 'TestProject'.

Default Execution queue is chosen for the database.

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

var InitalValueTS = db.SaveTimestep(VelX); // further fields an be 
//                                                  appended

BoSSSshell.WorkflowMgm.Sessions

#0: TestProject	InitialValueSession	12/05/2023 01:53:30	8d9a2627...

/// Now, we can use this timestep as a restart-value for the simulation:

C.SetRestart(InitalValueTS);