MPh

Pythonic scripting interface for Comsol Multiphysics®

Comsol is a commercial software application that is widely used in science and industry for research and development. It excels at modeling almost any (multi-)physics problem by solving the governing set of partial differential equations via the finite-element method. It comes with a modern graphical user interface to set up simulation models and can be scripted from Matlab® or via its native Java API.

This library brings the dearly missing power of Python to the world of Comsol. It leverages the universal Python-to-Java bridge provided by JPype to access the native API, and wraps it in a layer of pythonic ease-of-use. The Python wrapper only covers common scripting tasks, such as loading a model from a file, modifying some parameters, running the simulation, to then evaluate the results. Though the full functionality is available to those willing to drop down to the Java layer underneath.

Comsol models are marked by their .mph file extension, which stands for multi-physics. Hence the name of this library. It is open-source and in no way affiliated with Comsol Inc., the company that develops and sells the simulation software.

Installation

MPh is available on PyPI and can be readily installed via pip:

pip install mph

Add --user to the above command to make it a per-user installation, instead of system-wide, which may or may not be preferable. Run

pip uninstall mph

in order to remove the library from your system.

Requires JPype for the bridge from Python to Comsol’s Java API and NumPy for returning (fast) numerical arrays. pip makes sure the two Python dependencies are installed and adds them if missing.

Comsol, obviously, you need to license and install yourself. Versions from Comsol 5.3 onward are expected to work. If you want to use the client–server mode, make sure to once run comsolmphserver from a console first to enter the obligatory user name and password.

Support for Linux is limited: stand-alone clients do not work out of the box. You have to, once and for all, configure your shell’s environment to work around this issue. MacOS support is completely untested. It may work, but probably won’t. Find more details in the chapter Limitations.

Note how a separate Java run-time environment is not required, as Comsol ships with one already built in.

Tutorial

To follow along with this tutorial in an interactive Python session, if you wish to do so, make sure you have downloaded the demonstration model “capacitor.mph” from this library’s GitHub repository, and saved it in the same folder from which you run Python.

It is a model of a non-ideal, inhomogeneous, parallel-plate capacitor, in that its electrodes are of finite extent, the edges are rounded to avoid excessive electric-field strengths, and two media of different dielectric permittivity fill the separate halves of the electrode gap. Running the model only requires a license for the core Comsol platform, but not for any add-on module beyond that.

Starting Comsol

In the beginning was the client. And the client was with Comsol. And the client was Comsol. So let there be a Comsol client.

>>> import mph
>>> client = mph.Client(cores=1)

It takes roughly ten seconds until the client is initialized. In this example, the Comsol back-end is instructed to use only one processor core. If the optional parameter is omitted, it will use all cores available on the machine. Restricting this resource is useful when running several simulations in parallel. Note, however, that due to limitations of this library’s underlying Python-to-Java bridge, you cannot instantiate the Client class more than once. If you wish to accomplish that, in order to realize the full compute power of your simulation hardware, you need to start separate Python sessions, one for each client.

Managing models

Now that we have the client up and running, we can tell it to load a model file.

>>> model = client.load('capacitor.mph')

It returns a model object, i.e. an instance of the Model class. We will learn what to do with it further down. For now, it was simply loaded into memory. We can list the names of all models the client currently manages.

>>> client.names()
['capacitor']

If we were to load more models, that list would be longer. Note that the above simply displays the names of the models. The actual model objects can be recalled as follows:

>>> client.models()
[<mph.model.Model object at 0x000002CAF6C0BE80>]

We will generally not need to bother with these lists, as we would rather hold on to the model reference we received from the client. But to free up memory, we could remove a specific model.

>>> client.remove(model)

Or we could remove all models at once — restart from a clean slate.

>>> client.clear()
>>> client.names()
[]

Inspecting models

Let’s have a look at the parameters defined in the model.

>>> for parameter in model.parameters():
...     print(parameter)
...
parameter(name='U', value='1[V]', description='applied voltage')
parameter(name='d', value='2[mm]', description='electrode spacing')
parameter(name='l', value='10[mm]', description='plate length')
parameter(name='w', value='2[mm]', description='plate width')

Or the materials for that matter.

>>> model.materials()
['medium 1', 'medium 2']

They will be used by these physics interfaces:

>>> model.physics()
['electrostatic', 'electric currents']

To solve the model, we will run these studies:

>>> model.studies()
['static', 'relaxation', 'sweep']

Notice something? All features are referred to by their names, also known as labels, such as medium 1. But not by their tags, such as mat1, which litter not just the Comsol programming interface, but, depending on display settings, its graphical user interface as well.

Tags are an implementation detail. An unnecessary annoyance to anyone who has ever scripted a Comsol model from either Matlab or Java. Unnecessary because names/labels are equally enforced to be unique, so tags are not needed for disambiguation. And annoying because we cannot freely change a tag. Say, we remove a feature, but then realize we need it after all, and thus recreate it. It may now have a different tag. And any code that references it has to adapted.

This is Python though. We hide implementation details as much as we can. Abstract them out. So refer to things in the model tree by what you name them in the model tree. If you remove a feature and then put it back in, just give it the same name, and nothing has changed. You may also set up different models to be automated by the same script. No problem, as long as your naming scheme is consistent throughout.

Modifying parameters

As we have learned from the list above, the model defines a parameter named d that denotes the electrode spacing. If we know a parameter’s name, we can access its value directly.

>>> model.parameter('d')
'2[mm]'

If we pass in not just the name, but also a value, that same method modifies it.

>>> model.parameter('d', '1[mm]')
>>> model.parameter('d')
'1[mm]'

This particular model’s only geometry sequence

>>> model.geometries()
['geometry']

is set up to depend on that very value. So it will effectively change the next time it is rebuilt. This will happen automatically once we solve the model. But we may also trigger the geometry rebuild right away.

>>> model.build()

Running simulations

To solve the model, we need to create a mesh. This would also be taken care of automatically, but let’s make sure this critical step passes without a hitch.

>>> model.mesh()

Now run the first study, the one set up to compute the electrostatic solution, i.e. the instantaneous and purely capacitive response to the applied voltage, before leakage currents have any time to set in.

>>> model.solve('static')

This modest simulation should not take longer than a few seconds. While we are at it, we may as well solve the remaining two studies, one time-dependent, the other a parameter sweep.

>>> model.solve('relaxation')
>>> model.solve('sweep')

They take a little longer, but not much. We could have solved all three studies at once, or rather, all of the studies defined in the model.

>>> model.solve()

Evaluating results

Let’s see what we found out and evaluate the electrostatic capacitance, i.e. at zero time or infinite frequency.

>>> model.evaluate('2*es.intWe/U^2', 'pF')
array(1.31947349)

All results are returned as NumPy arrays. Though scalars such as this one could be readily cast to a (regular Python) float.

We could also ask where the electric field is strongest.

>>> (x, y, E) = model.evaluate(['x', 'y', 'es.normE'])
>>> E.max()
1455.5010798501924
>>> imax = E.argmax()
>>> x[imax], y[imax]
(-0.0005161744830262844, 0.004178390490229011)

Note how this time we did not specify any units. When left out, values are returned in default units. Here specifically, the maximum field strength in V/m and its coordinates in meters.

We also did not specify the dataset, even though there are three different studies that have separate solutions and datasets associated along with them. When not named specifically, the default dataset is used. That generally refers to the study defined first, here “static”. The default dataset is the one resulting from that study, here — inconsistently — named “electrostatic”.

>>> model.datasets()
['electrostatic', 'time-dependent', 'parametric sweep']

Now let’s look at the time dependence. The two media in this model have a small, but finite conductivity, leading to leakage currents in the long run. As the two conductivities also differ in value, charges will accumulate at the interface between the media. This interface charge leads to a gradual relaxation of the total capacitance over time. We can tell that from its value at the first and last time step.

>>> C = '2*ec.intWe/U^2'
>>> model.evaluate(C, 'pF', 'time-dependent', 'first')
array(1.31947349)
>>> model.evaluate(C, 'pF', 'time-dependent', 'last')
array(1.48410629)

The “first” and “last” time step defined in that study are 0 and 1 second, respectively.

>>> (indices, values) = model.inner('time-dependent')
>>> values[0]
0.0
>>> values[-1]
1.0

Obviously, the capacitance also varies if we change the distance between the electrodes. In the model, a parameter sweep was used to study that. These “outer” solutions, just like the time-dependent “inner” solutions, are referenced by indices, i.e. integer numbers, each of which corresponds to a particular parameter value.

>>> (indices, values) = model.outer('parametric sweep')
>>> indices
array([1, 2, 3], dtype=int32)
>>> values
array([1., 2., 3.]
>>> model.evaluate(C, 'pF', 'parametric sweep', 'first', 1)
array(1.31947349)
>>> model.evaluate(C, 'pF', 'parametric sweep', 'first', 2)
array(0.73678503)
>>> model.evaluate(C, 'pF', 'parametric sweep', 'first', 3)
array(0.52865545)

Then again, with a scripting interface such as this one, we may as well run the time-dependent study a number of times and change the parameter value from one run to the next. General parameter sweeps can get quite complicated in terms of how they map to indices as soon as combinations of parameters are allowed. Support for this may therefore be dropped in a future release — while the API is still considered unstable, which it is for as long as the version number of this library does not start with a 1 —, just to keep things simple and clean.

Saving results

To save the model we just solved, along with its solution, just do:

>>> model.save()

This would overwrite the existing file we loaded the model from. To avoid this, we could specify a different file name.

>>> model.save('capacitor_solved')

The .mph extension will be added automatically if it is not included in the first place.

Maybe we don’t actually need to keep the solution and mesh data around. The model was quick enough to solve, and we do like free disk space. We would just like to be able to look up modeling details somewhere down the line. Comsol also keeps track of the modeling history: a log of which features were created, deleted, modified, and in which order. Typically, such details are irrelevant. We can prune them by resetting that record.

>>> model.clear()
>>> model.reset()
>>> model.save('capacitor_compacted')

Most functionality that the library offers is covered in this tutorial. The few things that were left out can be gleaned from the API documentation.

Limitations

Java bridge

This library depends on the Python-to-Java bridge JPype. It therefore inherits its limitation, in that only one Java virtual machine can be managed within the same Python session, and thus only one Comsol client. If several simulations are to be run in parallel, distributed over independent processor cores in an effort to achieve maximum speed-up of a parameter sweep, they have to be started as separate Python (sub-)processes. This is a feasible work-around, but a limitation nonetheless.

For the same reason, the unit tests that come with MPh fail when collected and run via the testing framework pyTest. They must be run directly from the command line. Since they cannot be used for continuous integration testing anyway, given that they depend on Comsol being installed on the local machine, this is but a minor inconvenience.

Furthermore, there are some known, but unresolved issues with JPype’s shutdown of the Java virtual machine. Most notably, pressing Ctrl+C to interrupt an ongoing operation will usually crash the entire Python session. So do not rely on that to work.

(There is an alternative Java bridge, pyJNIus, which is not limited to one virtual machine, but then fails in another regard: A number of Java methods exposed by Comsol are inexplicably missing from the Python encapsulation.)

Creating models

The API is intentionally scope-limited to automating the simulation workflow, such as running parameter sweeps or optimization routines with customized, Python-powered post-processing. Exposing any and all Comsol features to create or alter every aspect of a model, and thus replicating the entire Java API, albeit in a more pythonic way, is out of scope.

Though users who do want to dig deeper may access the “pythonized” Java layer directly, via the .java attribute of Client instances (mapping to Comsol’s ModelUtil) as well as Model (mapping to Comsol’s model). Refer to Comsol’s programming manual for details.

This might even be worthwhile for Java developers of Comsol applications, as the interactive Python prompt provides easy introspection of the object hierarchy. In these circumstances, one might find the convenience function mph.inspect helpful, as it not only lists an object’s methods in a more readable way than Python’s built-in dir, but also displays the “properties” possibly defined on a model node.

In that event, it may equally help to call the returned Python model object something like pymodel and assign the name model to pymodel.java. Then you can just copy-and-paste Java (or even Matlab) code from the Comsol Programming Manual or as exported from the Comsol front-end. Python will gracefully overlook gratuitous semicolons at the end of statements, so this approach would even work for entire blocks of code. Keep in mind, however, that JPype cannot perform all type conversions silently in the background. Occasionally, when there is ambiguity in overloaded methods, you will have to cast types explicitly. Refer to the JPype documentation for help.

Client-server mode

In addition to running a stand-alone Comsol client, MPh also supports client–server connections. Client and server may either run on the same machine (which may or may not have benefits in terms of memory overhead) or they may communicate across the network. Note, however, that the latter scenario has never actually been tested — for lack of a network license required for that purpose. It may work as is, or may require transmitting login details: user name and password. File a GitHub issue if you are in a position to test this and want the problem solved.

Linux support

On Linux, the stand-alone Comsol client (as opposed to a client connecting to a server) will fail to start if all you did was install MPh. You will also have to add the names of shared-library folders to the environment variable LD_LIBRARY_PATH and explicitly export it. Otherwise initialization of the stand-alone client will throw a java.lang.UnsatisfiedLinkError because required external libraries cannot be found.

Add the following lines at the end of your shell configuration file, such as .bashrc:

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:\
/usr/local/comsol56/multiphysics/lib/glnxa64:\
/usr/local/comsol56/multiphysics/ext/graphicsmagick/glnxa64

Clearly, these paths depend on the installed Comsol version and will have to be adapted accordingly.

Requiring this variable to be set correctly limits the possibility of selecting a specific back-end from within MPh, as adding multiple installations to that search path would very likely lead to name collisions. Ideally, starting a stand-alone client would work out of the box, without the above manual configuration. However, this issue is currently unresolved.

macOS support

Support for macOS is experimental and entirely untested. The implementation is nearly identical to the one for Linux, and only deviates wherever the Comsol documentation indicates a difference in naming. As such, the environment variable for the library search path is instead called DYLD_LIBRARY_PATH, and Comsol would typically be found in, e.g., /Applications/COMSOL56/Multiphysics. Adjust your shell configuration accordingly. And please report back on the corresponding Github issue if you were able to test this.

API

Code documentation of the application programming interface provided by this library.

Model	Represents a Comsol model.
Client	Manages the Comsol client instance.
Server	Manages a Comsol server instance.
tools	Tools for software development.
discovery	Discovers Comsol installations on the local machine.