The 'v', 'x', and 'positiveVelocity' in spatialDistribution

[henrikt-ma]

Understanding the design of spatialDistribution isn't easy, see https://specification.modelica.org/master/operators-and-expressions.html#spatialdistribution. Here, I'm bringing up one small detail and one question.

Small detail: In the following, the v is not the name of any variable in the calling syntax:

where in0, in1, out0, out1, x, v are all subtypes of Real

Question: What is the rationale behind requiring the positiveVelocity to be provided separately; why can't the function just look at the change in x?

Edit: One more question: Should the resolution of initialPoints be considered the standardized way of controlling the internal resolution used in implementations, or should tools provide other means for users to control the internal resolution?

[HansOlsson]

Understanding the design of spatialDistribution isn't easy, see https://specification.modelica.org/master/operators-and-expressions.html#spatialdistribution. Here, I'm bringing up one small detail and one question.

Small detail: In the following, the v is not the name of any variable in the calling syntax:

where in0, in1, out0, out1, x, v are all subtypes of Real

Question: What is the rationale behind requiring the positiveVelocity to be provided separately; why can't the function just look at the change in x?

The rationale is that you sometimes want to consider an one-directional flow; so even if some numerical noise makes v<=0 you still want to view it as if v>0. A similar case is mentioned in the non-normative text.

Note that v is used in the text, even though it is not part of the calling syntax.

Edit: One more question: Should the resolution of initialPoints be considered the standardized way of controlling the internal resolution used in implementations, or should tools provide other means for users to control the internal resolution?

For the initial points it makes sense to use that (you cannot have higher resolution, and lower resolution seems bad), but for later values it is more complicated.

[henrikt-ma]

Edit: One more question: Should the resolution of initialPoints be considered the standardized way of controlling the internal resolution used in implementations, or should tools provide other means for users to control the internal resolution?

For the initial points it makes sense to use that (you cannot have higher resolution, and lower resolution seems bad), but for later values it is more complicated.

This is also what I would have guessed. I would then suggest that some text is added to clarify that the initialPoints resolution is independent of the resolution used in tools for approximation of the infinite-dimensional state.

[henrikt-ma]

Question: What is the rationale behind requiring the positiveVelocity to be provided separately; why can't the function just look at the change in x?

The rationale is that you sometimes want to consider an one-directional flow; so even if some numerical noise makes v<=0 you still want to view it as if v>0. A similar case is mentioned in the non-normative text.

Then it seems that the effect of inconsistent specification of (the change in) x and positiveVelocity deserves more clear normative specification, rather than just relying on a fairly bit piece of pseudo-code to speak for itself. I would have preferred something of more mathematical nature over the pseudo-code.

[henrikt-ma]

Note that v is used in the text, even though it is not part of the calling syntax.

Yes, I noted that, but didn't find it a valid reason to present it as being part of the calling syntax.

[casella]

Note that v is used in the text, even though it is not part of the calling syntax.

Yes, I noted that, but didn't find it a valid reason to present it as being part of the calling syntax.

But it isn't. The rationale of this operator is to provide native support for the solution of the advection equation, which is used for plug flow modelling and for wave propagation modelling, using the method of lines. So, the specification starts by recalling the advection equation, that uses v in its standard definition. It is then stated that v = der(x), and only then the calling syntax is given, which uses x.

I can't see anything wrong in this presentation.

[casella]

Question: What is the rationale behind requiring the positiveVelocity to be provided separately; why can't the function just look at the change in x?

Because that would go against giving a declarative definition of the operator, which I understand is what you are asking for as "something of more mathematical nature over pseudo code". The mathematical definition is stated up front: the operator describes the solution of a PDE, and has as inputs everything you need to compute that. How you compute it is up to the implementor.

Now, there is no such thing as "the change of x" in a Modelica operator, for the same reason that there is no operator in Modelica that can explicitly access the value at the previous integration step of a continuous-time variable to do stuff. The infamous "memory" operator of Simulink. No way. Allowing this would break the declarative nature of the language and open Pandora's box, the gates of hell, and so on 😄

Since the velocity v is not an input of the operator (it's none of its business to solve v = der(x)!) but you need to know its sign to determine the evolution of the solution, then the sign of v (and only the sign) must be an input of the operator. This is a mathematical formulation. How you actually implement it is a different story, and BTW the pseudo-code implementation is by definition not mandatory, because it is not formally defined.

[casella]

@AnHeuermann, you recently implemented the spatialDistribution() operator in OpenModelica, maybe you want to comment on that as well.

[eshmoylova]

My problem is with x and y used there. In one instance in the non-normative section at the very beginning it says z(x,t) everywhere else it is z(y,t). The PDE formula uses z(y,t) and v(t) v is given as derivative of x but x is not defined. Yet, the parameter specification for spatialDistibution makes use of x and not y. Whenever z is mentioned, it has y and when spatialDistribution is defined it is x. Are y and x supposed to be the same? If not, what is x and why isn't y present in the parameters of the spatialDistribution.

[casella]

@eshmoylova you are absolutely right, that is of course a mistake, we should have either x or y throughout the text. In fact, this issue has been dangling there for a while, see e.g. #1729.

I dug out some documents. The discussion on this operator started in 2009, based on the transportDelay operator of Dymola, see #172. This document was later discussed in the Fluid Group at the 69th Modelica Design Meeting in Munich, Jan 2011 Minutes-Fluid-Group.pdf, and later finalized at the 72nd Modelica Design Meeting in Linkoping in Sep 2011 spatialDistribution.pdf. In both cases the text used x throughout. For some reason, however, the text that got into Modelica 3.3 had y instead of x.

Based on the comments given in #1729, I guess at some points we had some documents/examples where we both had a spatial coordinate x (with unit "m") and a non-dimensional normalized spatial coordinate y = x/L. We eventually opted for the choice of the normalized spatial coordinate, to avoid further complications with units, and to avoid any reference to the non-normalized spatial variable, for conciseness. I guess the use of y in parts of the specification is a leftover of that process, even though the Design Meeting Minutes always used x.

Bottom line, it's finally time to put x everywhere, as in the original minutes.

[casella]

See PR #3113

[AnHeuermann]

@AnHeuermann, you recently implemented the spatialDistribution() operator in OpenModelica, maybe you want to comment on that as well.

I think myself and @kabdelhak did find the standard on this confusing as well. The pseudo-code part was helpful, but at the same time confusing as well. Problem was that I pictured a conveyor belt in my head, but that doesn't capture everything from spatialDistribution.

[casella]

I think myself and @kabdelhak did find the standard on this confusing as well. The pseudo-code part was helpful, but at the same time confusing as well. Problem was that I pictured a conveyor belt in my head, but that doesn't capture everything from spatialDistribution.

Could you be more specific on the confusing parts, and maybe suggest how to improve it?

As to the conveyor belt image (or plug flow in a pipe), of course it captures the mathematical problem, but not how its solution can be programmed in software.

[AnHeuermann]

We did this around a year ago. Sorry, I don't remember anything any more.

But I think the concept of

  Real points[:];
  Real values[:];

being storage over multiple function calls is strange and doesn't transalte well to Modelica.

And we discussed about initialization: What are allowed values for initPoints?

In the end we sayed something like {0.0, 0.1, 0.2, 0.2, 0.3, 0.6, 1.0} would mean that there is an event stored at 0.2.

Also we implemented it in two different functions:

1. storeSpatialDistribution --> Store something new on the spatial distribution
2. spatialDistribution --> Remove something (output) from the spatial distribution on an accepted step

Because OpenModelica is allowed to fail at any point and its easy to undo storeSpatialDistribution, but not the removal of stuff from spatial distribution.

So from the pseudo-code example I would remove points and values or switch to a language (C maybe) where you can have this as internal member of the spatialDistribution operator.

[HansOlsson]

So from the pseudo-code example I would remove points and values or switch to a language (C maybe) where you can have this as internal member of the spatialDistribution operator.

Pseudo-code without points and values would be too limited and not fully explain it. Writing C-code might be good, if someone volunteered - and there are no obvious bugs. It might also be that another language (C++ might be an obvious choice) would be better by avoiding the low-level memory-handling of arrays.

[HansOlsson]

Moving it to a future mile-stone as it is not urgent, and someone needs to write it.

[casella]

@HansOlsson, we have an open-source implementation of spatialDistribution in OpenModelica, see spatialdistribution.h and spatialdistribution.c, written by @AnHeuermann and @kabdelhak, that we could contribute to the specification as a reference implementation.

It is quite complete, as it includes an adaptive resizable circular buffer and propagation of events from inputs to outputs. It has been successfully validated against non-trivial test cases with known analytic closed-form solution.

However, this seems a bit too involved for me to find place in the specification. Maybe we could add a link to it in the specification text?

[HansOlsson]

However, this seems a bit too involved for me to find place in the specification. Maybe we could add a link to it in the specification text?

Adding such implementation hints seems a bit outside the scope of the specification; and we might soon be over-burdened.

[AnHeuermann]

Adding such implementation hints seems a bit outside the scope of the specification; and we might soon be over-burdened.

I agree and the implementation is using OpenModelica specific data types and structures. And with nearly 1000 LOC it's way to long to act as an easy to understand reference.

[casella]

I also agree.

@AnHeuermann, what you could do instead is to add the test case to ModelicaTest. That would be useful for all tool vendors to check the correctness of their results, since the model also computes the exact closed-form solution.

[HansOlsson]

It seems it can now be closed (and there was a clarification in #3113 ). The test-suite for the semantics is https://github.com/modelica/Modelica-Compliance - generally they have been more focused on semantic restrictions, but it makes perfect sense to also add a check for this.