CESDSOL

Complex Equations — Simple Domain partial differential equations SOLver

CESDSOL is a C++ library providing capabilities to solve arbitrary partial differential equations with emphasis on highly nonlinear problems arising in field theory. It is battle-tested on numerous problems of finding stationary configurations for field-theoretic and condensed matter solitons, including gravitating ones, and black holes interacting with matter fields and of studying dynamics of soliton collisions.

Features

• Stationary and explicit transient problems
• Arbitrary spatial dimension of problems
• Arbitrary number of equations
• Arbitrary order of derivatives, both spatial and temporal
• Arbitrarily-spaced cartesian grids
• Arbitrary boundary conditions
• Finite-difference discretization of arbitrary order of accuracy
• Field equations can be supported with discrete variable equations
• Variety of linear, nonlinear and ODE solvers
• Intel MKL as linear algebra backend
• Provided with Wolfram Language interface for generation of problem descriptor and reading of solver output

CESDSOL is developed in a very modular fashion allowing for extensions in different directions. Degrees of freedom include:

• New problem types
• New grid types
• New discretization types
• New math backends
• New solvers

Building

CESDSOL is just an almost header-only library so it's up to user to integrate it into application, be it static or dynamic library or executable. CESDSOL is provided with cmake file for typical mode of operation where for each problem new executable is built. In particular it can be used to build examples. For example, use command cmake --build . -DMainPath=examples/StaticSineGordon/StaticSineGordon.cpp in order to build static Sine-Gordon kink example. Generally -DMainPath is used to provide .cpp file with main() function where problem definition and solution process is set up.

CESDSOL is written in C++20, so C++20 compliant compiler is required. It's tested only on Windows, so Visual Studio of version not less than 16.10 is recommended. Currently the only supported math backend is Intel MKL, so it's also required for the best experience.

Usage

Problem setting in CESDSOL is a pretty complex and tedious task, so Mathematica interface is recommended. In order to use CESDSOL interface from Mathematica first you have to load packgage CESDSOL.wl from Mathematica folder. Three routines are used for work with CESDSOL.

CESDSOL`StationaryProblem[eqs, vars, coords, grid, options]

is used to generate problem input for stationary problems, where eqs are left-hand sides of equations (in form f(x)=0) including boundary conditions, vars are fields, coords are spatial coordinates and grid is a grid to solve problem on. Options include the following:

• CalculateJacobian — if jacobian calculation is needed. True by default, can be switched to False if for problem solution methods not requiring jacobian will be used
• Parameters — list of problem parameters
• DiscreteEquations — list of discrete variable equations
• Variables — list of discrete variables
• Name — problem name
• LocalPIEs — local parameter-independent expressions. Can be used to alias some repeatedly used expressions which depend only on coordinates
• GlobalPIEs — global parameter-independent expressions. Can be used to alias some repeatedly used expressions which do not depend on coordinates
• LocalVIEs — local variable-independent expressions. Can be used to alias some repeatedly used expressions which depend on coordinates and parameters
• GlobalVIEs — global variable-independent expressions. Can be used to alias some repeatedly used expressions which depend only on parameters
• LocalVDEs — local variable-dependent expressions. Can be used to alias some repeatedly used expressions which can depend on fields
• GlobalVDEs — global variable-dependent expressions. Can be used to alias some repeatedly used expressions which can depend on dicrete variables
• Integrals — integrals appearing in equations
• Discretization — discretization to use. Currently only StructureFiniteDifferenceDiscretization[order] is used, where order is accuracy order finite-difference of discretization
• LocalOutput — local output expressions
• GlobalOutput — global output expressions
• IntegralOutput — integral output expressions
• PointOutput — point-interpolated output expressions

For example for static Sine-Gordon problem problem definition is generated by the following call:

CESDSOL`StationaryProblem[{
  {
   phi''[x] - g Sin[phi[x]],
   Left[x] -> phi[x] - π,
   Right[x] -> phi[x] - 2 π
   }
  },
 {phi},
 {x},
 DirectProductGrid[UniformRange[0, 10, 1000]],
 Parameters -> {g},
 Name -> "StaticSineGordon"]

Note that equations are provided as list

{
 {
  phi''[x] - g Sin[phi[x]],
  Left[x] -> phi[x] - π,
  Right[x] -> phi[x] - 2 π
  }
 }

where each element is list including left-hand side of equation phi''[x] - g Sin[phi[x]] and its boundary condition on left side of x-range Left[x] -> phi[x] - π and boundary condition on right side of x-range Right[x] -> phi[x] - 2 π.

CESDSOL`ExplicitTransientProblem[eqs, vars, coords, grid, options]

is used to generate problem input for epxlicit transient problems, where eqs are equations in form highest fild time derivative -> {right-hand sides of equations}. It supports all the options of StationaryProblem except CalculateJacobian. As an example problem definition for fermion-kink collision is generated by the following call:

ExplicitTransientProblem[{
  Derivative[0, 2][f][x, t] -> {Derivative[2, 0][f][x, t] - Sin[f[x, t]] + 2 g (u[x, t] v[x, t] + p[x, t] q[x, t])}, 
  Derivative[0, 1][u][x, t] -> {-Derivative[1, 0][q][x, t] + g f[x, t] q[x, t]},
  Derivative[0, 1][v][x, t] -> {Derivative[1, 0][p][x, t] + g f[x, t] p[x, t]},
  Derivative[0, 1][p][x, t] -> {Derivative[1, 0][v][x, t] - g f[x, t] v[x, t]},
  Derivative[0, 1][q][x, t] -> {-Derivative[1, 0][u][x, t] - g f[x, t] u[x, t]}
  },
 t, {x},
 DirectProductGrid[UniformRange[-100, 100, 2000]],
 Parameters -> {g},
 LocalOutput -> {"norm" -> Sqrt[u[x, t]^2 + v[x, t]^2 + p[x, t]^2 + q[x, t]^2]},
 Name -> "FermionSineGordonKinkCollision"]

When problem definition is generated just paste it into C++ code, set up solver and run the code! There is a variety of solvers. The recommended setup for stationary problems is ModifiedNewton solver with MKL::PARDISO as linear solver and GoldenSectionSearch as line searcher:

auto pardiso = std::make_unique<MKL::PARDISO<double>>();
auto gss = MakeLineSearcher<GoldenSectionSearch>(*problem);
auto newton = MakeNonlinearSolver<ModifiedNewton>(*problem, std::move(pardiso), std::move(gss));

Use then newton.Solve(*problem) to solve the problem for current values of parameters or use FixedStepParametericSweeper of AdaptiveParametericSweeper to automatically sweep over parameters.

The final step is data extraction. User can of course set up any suitable serialization, but default Seializer is recommended. It generates data in binary format, including grid, parameters, variable and output values. Mathematica interface includes routine

CESDSOL`ReadSolution[path]

for data import for further analysis.

Refer to examples from examples folder for further information.