PLaSK library

Loading...
Searching...
No Matches
Solvers

# About solvers

Solvers in PLaSK are calculations units. Each is represented by subclass of plask::Solver.

Each solver has one or more input and one or more output realized by providers and receivers mechanism. This allows for communication between solvers.

Typically, each solver is connected with geometry object which describe physical property (mainly materials) of laser or its fragment. This geometry object defines (local) calculation space for solver.

Typically, each solver also includes mesh which represent set of point (in calculation space) in which solver calculate its result. If another solver requests for data in points other than these includes in mesh, result can be interpolated.

# How solvers are typically used?

TODO

Note that typically, solvers are used from python scripts.

# How to write a new calculation solver?

To write solver you should:

1. If you write an official solver for PLaSK, create a subdirectory structure in directory solvers. The subdirectories in solvers should have a form solvers/category/your_solver, e.g. solvers/optical/FDTD. It is also possible to write fully external solvers by PLaSK users, using the descpription below (TODO).
2. Copy solvers/skel/CMakeLists.txt from to your subdirectory. You may edit the copied file to suit your needs.
3. By default all the sources should be placed flat in your solver subdirectory and the Python interface in the subdirectory your_solver/python.

Once you have your source tree set up, do the following:

1. Write new class which inherits from plask::SolverOver<SPACE_TYPE> or plask::SolverWithMesh<SPACE_TYPE, MESH_TYPE>. SPACE_TYPE should be one of Geometry2DCartesian, Geometry2DCylindrical, or Geometry3D and should indicate in what space your solver is doing calculations. If you want to allow user to specify a mesh for your solver, inherit from plask::SolverWithMesh<SPACE_TYPE, MESH_TYPE>, where MESH_TYPE is the type of your mesh.
2. Implement plask::Solver::getClassName method. This method should just return the pretty name of your solver class (the same you use in Python and XML).
3. Implement plask::Solver::onInitialize and optionally plask::Solver::onInvalidate methods.
4. Place in your class public providers and receivers fields:
   - provider fields allow your solver to provide results to another solvers or to reports,
   - receiver fields are used to getting data which are required for calculations from another solvers (precisely, from its providers),
   - you don't have to care about connection between providers and corresponding receivers (this is done externally),
   - note that most providers are classes obtained by using plask::ProviderFor template,
   - more details can be found in Providers and receivers.
5. If you need boundary conditions, place in your class public plask::BoundaryConditions fields which are containers of boundary-condition pairs.
6. Implement loadConfiguration method, which loads configuration of your solver from XML reader.
7. Typically, implement calculation method. This method is a place for your calculation code. You don't have to implement it if you don't need to do any calculations. You can also write more methods performing different calculations, however, you need to clearly document them. Each calculation method must call plask::Solver::initCalculation() as the first operation.
8. Optionally implement plask::Solver::getClassDescription method. This method should just return description of your solver.
9. Write the Python interface to your class using Boost Python. See the Boos Python documentation or take a look into `solvers/skel/python/solver.cpp` (for your convenience we have provided some macros that will facilitate creation of Python interface).
10. Write `solvers.yml` file that documents `xpl` configuration of your solver for generating the `xpl` reference in the user manual and for automatic configuration panel creation in GUI.
11. Finally write a good user manual for your solver ;)

## Writing solvers in depth

Below we explain the above steps in detail on a simple example. When following this tutorial work already on your own solver, or write a sample code in a separate directory in order to keep the PLaSK source tree (and repository) clean.

Assume that we want to write a solver computing a waveguide effective index and optical field intensity for edge emitting lasers. The solver performs its computation using finite differences method. Hence, we name it FiniteDifferencesSolver.

To start we create a subdirectory solvers/optical/finite_diff under the PLaSK trunk write a file finite_differences.h to it. Also we need to copy the CMakeLists.txt file from solvers/skel to our directory. In most cases we will only need to edit the line with the command `project` in this file, unless our solvers uses some external libraries. If for example, our solver uses LAPACK, the CMakeLists.txt should look like this (In the below example all comments are skipped):

The project name should math the pattern plask/solvergroup/solverlib, so in our case it will look like plask/optical/finite_diff.

    project(plask/optical/finite_diff)

    set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../../../cmake)
    find_package(PLaSK)

    find_package(LAPACK)
    set(SOLVER_LINK_LIBRARIES ${LAPACK_LIBRARIES})

    make_default()

If you are working within the PLaSK source, remember to add all the directories and files to the subversion repository (only for your real solver naturally)!

### Solver C++ class

Now we assume that our solver uses rectilinear mesh provided by the user, so we inherit our class from plask::SolverWithMesh<plask::Geometry2DCartesian, plask::RectangularMesh<2>>. So our header finite_differences.h should begin as follows:

#ifndef PLASK__SOLVERS__OPTICAL__FINITE_DIFFERENCES_H
#define PLASK__SOLVERS__OPTICAL__FINITE_DIFFERENCES_H
#include <plask/plask.hpp>
namespace plask { namespace solvers { namespace optical_finite_differences { // put everything in private namespace
class PLASK_SOLVER_API FiniteDifferencesSolver: public plask::SolverWithMesh < plask::Geometry2DCartesian, plask::RectangularMesh<2> >
{

(note that PLASK_SOLVER_API is required for proper exporting/importing the solver to/from DLL file)

Then, you declare all the fields and methods of the class. We will skip all the fields that are required privately for computations from this tutorial and focus only on the ones necessary for PLaSK interface. Assume, that FiniteDifferencesSolver needs a temperature distribution and wavelength as an input, and outputs effective index of the waveguide and the optical field intensity. Additionally, boundary conditions of the first kind on temperature is needed. Hence, we declare the following providers, receivers and boundary conditions:

public:
plask::ReceiverFor<plask::Wavelength> inWavelength;
plask::ReceiverFor<plask::Temperature, plask::Geometry2DCartesian> inTemperature;
plask::ProviderFor<plask::EffectiveIndex>::WithValue outNeff;
plask::ProviderFor<plask::LightMagnitude, plask::Geometry2DCartesian>::Delegate outLightMagnitude;
plask::BoundaryConditions<plask::RectilinearMesh2D::Boundary, double> boundaryConditionsOnField;

In the code above, we have declared two receivers (by convention in PLaSK, names of every receiver in all solvers should begin with in prefix). `inWavelength` can receive a single value either from some connected solver (e.g. a computed gain maximum) or specified by the user. On the other hand `inTemperature` receives temperature distribution, which is a scalar field. For this reason it is necessary to specify the type of the space in which this distribution will be computed. Naturally it should match the working space of your solver.

When declaring providers, one needs to specify how the provided value can be obtained. In our case `outNeff` (again, the names of every provider in all solvers should begin with out prefix) will hold its value internally i.e. you will be able to assign a value to it as simply as `outNeff = 3.5;`. On the other hand, `outLightMagnitude` is a delegate provider i.e. you will need to write a method which computes the light intensity on demand (we will later show you how).

As your solver inherits plask::SolverWithMesh, there is already a shared pointer to the `geometry` and `mesh` available. However, it might be a good idea to create the class member field for the light intensity computed on the solver mesh to be able to efficiently provide it to other solvers. You can use any type for array for this (e.g. std::vector<double>), but—as the solvers exchange field data using plask::DataVector class—it is best to use this class for internal storage as well (it behaves mostly like std::vector, with some additional improvements required for shared memory management between the solvers). So we declare the private field:

private:
plask::DataVector<double> computed_light_intensity;

Now, you can write the constructor to your class. By convention this constructor should take no configuration arguments as all the solver configuration parameters must be able to be set by the user afterwards (in future you might want to create a configuration window for your solver for GUI). The only constructor parameter is the name, which can be provided by user, and which should be passed to the parent class. In addition, you should write getClassName method, which returns the category and the name of your solver class as seen by the end user (it does not need to be the same as your real class name, but it should match the Python class name and solver name in XML file).

public:
FiniteDifferencesSolver(const std::string name& name=""):
plask::SolverWithMesh<plask::Geometry2DCartesian,plask::RectilinearMesh2D>(name),
outLightMagnitude(this, &FiniteDifferencesSolver::getIntensity) // attach the method returning the light intensity to delegate provider
{
inTemperature = 300.; // set default value for input temperature
}
virtual std::string getClassName() const { return "optical.FiniteDifferences2D"; }

In the above illustration, we initialize the `outLightMagnitude` provider with the pointer to the solver itself and the address of the method computing the light intensity (we write this method later). Also, we set the default value of the temperature to 300 K in the whole structure. As there is not default value for inWavelength, the user will have to provide it (either manually or from some wavelength provider) or the exception will be raised when we try to retrieve the wavelength value in our computation method.

Before we write a computation method (or several computation methods), we must realize that the solver can be in two states: initialized or invalidated. When the user creates the solver object, sets the geometry and the mesh, there is still no memory allocated for our results (i.e. `computed_light_intensity` is an empty vector) nor the results are know. This is called an invalidated state. When in this state, the first step before any calculation is to allocate the memory: resize the `computed_light_intensity` to the proper size, allocate some internal matrices (not mentioned in this tutorial), etc. The size of the allocated memory depends e.g. on the mesh size, so it cannot be done in the constructor. For this reason, in the beginning of every calculation function you should call the method initCalculation(), which is defined in plask::Solver base class. If the solver has been in invalidated state, this it will call virtual method onInitialize(), in which you can put all your initialization code. Then the solver will be put in initialized state (hence, subsequent calls to computational method will not call `onInitialize()` unless solver is forcefully invalidated (see below). The code on the initialization method may look like this:

protected:
virtual void onInitialize() {
if (!geometry) throw NoGeometryException(getId()); // test if the user has provided geometry
if (!mesh) throw NoMeshException(getId()); // test if the user has provided the mesh
// do any memory allocations or your internal matrices, etc.
}

Even after some computations have been performed, the user might change the geometry of the structure or the mesh. In such case the results of your computations becomes outdated and the sizes of some matrices might need to change. In such situation, the solver is put back to invalidated state and, similarly to the initialization, the virtual method onInvalidate() is called. This method may look like this:

virtual void onInvalidate() {
outNeff.invalidate(); // clear the computed value
computed_light_intensity.reset(); // free the light intensity array
// you may free the memory allocated in onInitialize unless you use some containers with automatic memory management
}

This method can be also forcefully called by the user issuing `your_solver.invalidate();` command. This might be done in order to free them memory for some other computations. For this reason you should free all large chunks of memory in `onInvalidate()`. However, If you store your data in an array with no built-in memory management (e.g. old-style C-arrays), you have to check if the memory has been previously allocated, as `onInvalidate()` might be called before `onInitialize()`. Furthermore you should call this method from your class destructor in such case.

After the solver has been invalidated, before the next computations `onInitialize()` will be called, so the new memory will be allocated as needed.

Now you can write your core computational methods. Just remember to call initCalculation() in the beginning. Furthermore, if the provided values change (and they probably do), you must call method fireChanged() for each of your providers, to notify the connected receivers about the change. Here is the sample code:

public:
void compute() {
initCalculation(); // calls onInitialize if necessary and puts solver in initialized state
computed_light_intensity.reset(); // clear the previously computed light intensity (if any)
DataVector<double> temperature = inTemperature(*mesh); // gets temperature computed by some other solver
double wavelenght = inWavelength(); // gets the wavelength
// [...] perform your calculations
outNeff = new_computed_effective_index;
outNeff.fireChanged(); // inform others that we have computed a new effective index
outLightMagnitude.fireChanged(); // we will also be providing a new light intensity
}

Assume that in the above sample computation method, we did not compute the light intensity. We will do it only if someone wants to know it. Of course the choice when to compute what depends strongly on your particular solver and you should focus on efficiency (i.e. compute only what is necessary) when making decisions on this matter.

The last thing to do is to write the method called by the delegate provider `outLightMagnitude` when someone wants to get the optical field intensity distribution. The arguments and return value of this method depend on the provider type. For interpolated fields they will look like in the following example:

protected:
DataVector<const double> getIntensity(const plask::MeshD<2>& destination_mesh, plask::InterpolationMethod interpolation_method=INTERPOLATION_DEFAULT) {
if (!outNeff.hasValue()) throw NoValue(LightMagnitude::NAME); // this is one possible indication that the solver is in invalidated state
if (computed_light_intensity.size() == 0) // we need to compute the light intensity
{
computed_light_intensity.reset(mesh.size()); // allocate space for the light intensity
// [...] compute the light intensity
}
// automatically interpolate your data to the requested mesh
return interpolate(*mesh, computed_light_intensity, destination_mesh,
getInterpolationMethod<INTERPOLATION_LINEAR>(interpolation_method), this->geometry);
}

The important objects of the above method are the first and the last lines. In the former one, we check if the computations have been performed and are up-to-date (remember, we have cleared the value of `outNeff` in `onInvalidate()`). Otherwise we throw an exception. In the last line we use plask::interpolate function to interpolate our data to the receiver mesh (which is provided as `destination_mesh` argument).

Passing this->geometry helps to automatically consider mirror and periodic boundaries, so the requested points will be wrapped into your computational domain correctly. And getInterpolationMethod changes INTERPOLATION_DEFAULT method to some real one.

Our solver can perform computations now. However, if it has any configuration to load, we can read it from XML file. To do this, we should reimplement `loadConfiguration` method. It reads the configuration from the current XML file using plask::XMLReader, by walking through the consecutive tags. It is important to call `parseStandardConfiguration` if you encounter unknown tag. Below you have an example:

void loadConfiguration(XMLReader& reader, Manager& manager) {
while (reader.requireTagOrEnd()) {
if (reader.getNodeName() == "newton") {
newton.tolx = reader.getAttribute<double>("tolx", newton.tolx);
newton.tolf = reader.getAttribute<double>("tolf", newton.tolf);
newton.maxstep = reader.getAttribute<double>("maxstep", newton.maxstep);
newton.method = reader.enumAttribute<MethodEnumType>("method") // MethodEnumType is some declared enum
.value("raphson", METHOD_RAPHSON) // METHOD_RAPHSON is one of the enum values
.value("secant", METHOD_SECANT, 3) // METHOD_SECANT can be selected with "secant" or its first 3 letters ("sec")
.get(newton.method);
reader.requireTagEnd();
} else if (reader.getNodeName() == "wavelength") {
auto wavelength = reader.getAttributere<double>("value");
if (wavelength) inWavelength.setValue(*wavelenght);
} else if (reader.getNodeName() == "boundary") {
manager.readBoundaryConditions(source, boundaryConditionsOnField);
reader.requireTagEnd();
} else
parseStandardConfiguration(reader, manager, "<geometry>, <mesh>, <newton>, or <wavelength>");
}
}

In the above example we assume that we have some local `struct` with parameters of Newton algorithm: `tolx`, `tolf`, and `maxstep`. They are read from the corresponding attributes of a <newton> tag, with default value equal to the current value of the corresponding parameter. Furthermore, user can optionally specify a wavelength, which we set as a specified input value of the inWavelength receiver (single-value receivers can be connected to providers, however they can also have assigned value as normal variables, although in any case to read their values you must remember about using parenthesis, e.g. `w = inWavelength()`).

The XML file to read by the above method can look as follows (although you should rather use XML attributes to set simple parameters, in order to make the XML file consistent for all the solvers).:

    <optical lib="finite_diff" solver="FiniteDifferences2D">
        <newton tolx="0.0001" tolf="1e-9" maxstep="500" />
        <wavelenght>1000</wavelength>
    </optical>

You can now finish your class definition:

}; // class FiniteDifferencesSolver
}}} // namespace plask::solvers::optical_finite_differences
#endif // PLASK__SOLVERS__OPTICAL_FINITE_DIFFERENCES_H

### Writing Python interface

Once you have written all the C++ code of your solver, you should export it to the Python interface. To do this, create a subdirectory named python in your solver tree (i.e. solvers/optical/finite_diff/python) and copy your_solver.cpp from solvers/skel/python there (changing its name e.g. to finite_differences_python.cpp).

Then edit the file to export all necessary methods and public fields. The contents of this file should look more less like this:

#include <cmath>
#include <plask/python.hpp>
using namespace plask;
using namespace plask::python;
#include "../finite_differences.hpp"
using namespace plask::solvers::optical_finite_differences; // use your solver pricate namespace
BOOST_PYTHON_MODULE(fd)
{
{CLASS(FiniteDifferencesSolver, "FiniteDifferences2D",
"Calculate optical modes and optical field distribution using the finite\n"
"differences method in Cartesian two-dimensional space.")
RECEIVER(inWavelength, "Wavelength of the light");
RECEIVER(inTemperature, "Temperature distribution in the structure");
PROVIDER(outNeff, "Effective index of the last computed mode");
PROVIDER(outLightMagnitude, "Light intensity of the last computed mode");
BOUNDARY_CONDITIONS(boundary, boundaryConditionsOnField, "Boundary conditions of the first kind (constant field)");
METHOD(compute, compute, "Perform the computations");
}
}

`BOOST_PYTHON_MODULE` macro takes the name of the package with your solver (without quotation marks). Your solver will be accessible to Python if the user imports this package as:

import plask.optical.fd

The arguments of the `CLASS` macro are your solver class, the name in which it will appear in Python, and short solver documentation (mind the braces outside of the `CLASS:` they are important if you want to put more than one solver in a single interface file, so they will appear in a single package).

Next you define your exported class member fields, properties (fake fields, which call your class methods on reading or assignment), methods, providers, and receivers (you have probably noticed that providers and receivers are just class member fields, but they need to be exported using separate macros, due to some additional logic necessary). Below, there is a complete list of macros exporting class objects and we believe it is self-explanatory:

METHOD(python_method_name, method_name, "Short documentation", "name_or_argument_1", arg("name_of_argument_2")=default_value_of_arg_2, ...);
RO_FIELD(field_name, "Short documentation"); // read-only field
RW_FIELD(field_name, "Short documentation"); // read-write field
RO_PROPERTY(python_property_name, get_method_name, "Short documentation"); // read-only property
RW_PROPERTY(python_property_name, get_method_name, set_method_name, "Short documentation"); // read-write property
RECEIVER(inReceiver, "Short documentation"); // receiver in the solver
PROVIDER(outProvider, "Short documentation"); // provider in the solver

When defining methods that take some arguments, you need to put their names after the documentation, so the user can set them by name, as for typical Python methods. You can also specify some default values of some arguments as shown above.

After you successfully compile the solver and the Python interface (just run `cmake` and `make`), the user should be able to execute the following Python code:

import numpy
import plask
import plask.optical.fd
solver = plask.optical.fd.FiniteDifferences2D()
solver.geometry = plask.geometry.Cartesian2D(plask.geometry.Rectangle(2, 1, "GaN"))
solver.mesh = plask.mesh.Rectangular2D(Regular(0,2,100), Regular(0,1,100))
solver.inTemperature = 280
solver.compute()
print(solver.outNeff())

### Writing automatic solver tests

TODO

## Writing configuration description

Once your solver is working, you should create the solvers.yml file in the main directory of your solver. In this file you should specify all the XML configuration tags and attributes that are read from an `xpl` file. An example of such file can be found in `solvers/skel/solvers.yml`. For the above example the file should look like:

- solver: YourSolver
lib: your_solver
category: skel
geometry: YourGeometry
mesh: YourMesh
tags:
- tag: configuration-tag
label: Configuration Tag
help: Documentation of your configuration tag
attrs:
- attr: tag-attr
label: Tag attribute
type: attribute type
unit: unit
help: Pull all the configuration attributes like this.
- attr: other-attr
label: Other attribute
type: choice
choices:
- First choice
- Second choice
- Third choice
help: Another configuration attribute.
- tag: another-configuration-tag
label: Different Configuration Tag
help: ''
attrs:
- attr: attribute
label: Some attribute
type: attribute type
unit: unit
help: Pull all the configuration attributes like this.
- bcond: boundary_condtion
label: Boundary Condition
- bcond: other_boundary_condtion
label: Another Boundary Condition
help: Additional description of boundary condition in help.
providers:
- outProvider: QuantityIfNameDifferent
receivers:
- inReceiver

This concludes our short tutorial. Now you can go on and write your own calculation solver. Good luck!