Issue with MultiBody models of ScalableTestSuite with the NF

[Francesco Casella]

The last remaining models of ScalableTestSuite which are not handled by the NF are the FlexibleBeamModelica and StringModelica models. The issue is quite intricate but I think I understood what happens.

For example, you can check ​ScalableTestSuite.Mechanical.Strings.ScaledExperiments.StringModelica_N_2, which reports this error:

ScalableTestSuite_noopt_ScalableTestSuite.Mechanical.Strings.ScaledExperiments.StringModelica_N_2_06inz.c:4331:38:
error: use of undeclared identifier '$Pbodybox$lB1$rB$Pbody$PR_start$PT'
  copy_real_array_data_mem(tmp2._T, &$Pbodybox$lB1$rB$Pbody$PR_start$PT);

The issue is with the parameter-binding equation

  final parameter Frames.Orientation R_start=
      Modelica.Mechanics.MultiBody.Frames.axesRotations(
        sequence_start,
        angles_start,
        zeros(3));

of the Modelica.Mechanics.MultiBody.Parts.Body model, which is moved as a whole by the new front end into the initial equation section as

  bodybox[1].body.R_start = Modelica.Mechanics.MultiBody.Frames.axesRotations(
  {1, 2, 3}, bodybox[1].body.angles_start, {0.0, 0.0, 0.0});

which the back-end cannot handle. The old-front end constant-evaluates the binding equation and expands it to its scalar components

  final parameter Real bodybox[1].body.R_start.T[1,1] = 1.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.T[1,2] = 0.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.T[1,3] = 0.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.T[2,1] = 0.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.T[2,2] = 1.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.T[2,3] = 0.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.T[3,1] = 0.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.T[3,2] = 0.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.T[3,3] = 1.0 "Transformation matrix from world frame to local frame";
  final parameter Real bodybox[1].body.R_start.w[1](quantity = "AngularVelocity", unit = "rad/s") = 0.0 "Absolute angular velocity of local frame, resolved in local frame";
  final parameter Real bodybox[1].body.R_start.w[2](quantity = "AngularVelocity", unit = "rad/s") = 0.0 "Absolute angular velocity of local frame, resolved in local frame";
  final parameter Real bodybox[1].body.R_start.w[3](quantity = "AngularVelocity", unit = "rad/s") = 0.0 "Absolute angular velocity of local frame, resolved in local frame";

which is trivial for the back-end to handle.

However, looking at the axesRotation() function call, the parameter R_start is final, sequence_start has Evaluate = true and zeros(3) is a literal constant, but angles_start doesn't have Evaluate = true, so I'm not sure why the old front-end constant-evaluates the right-hand-side of the binding equation. Maybe just the fact of being final was used to trigger that, though it needs not necessarily be the case: a final parameter may depend on another parameter that is not evaluated and that should be changeable at runtime.

axesRotation has Inline = true, but I understand this is only relevant for the back-end, the front-end never inlines anything.

So, I don't think we should try to reproduce the old front-end behaviour here, because it performs an unnecessary evaluation. On the other hand, we need a viable solution to get ScalableTestSuite to run ASAP.

I checked ​Modelica.Mechanics.MultiBody.Examples.Elementary.DoublePendulum, which runs fine with the NF, and produces the initial equation

  boxBody1.body.R_start = Modelica.Mechanics.MultiBody.Frames.axesRotations(
  {1, 2, 3}, boxBody1.body.angles_start, {0.0, 0.0, 0.0});

which is handled correctly. So, I understand there is some problem with code generation of

  bodybox[1].body.R_start = Modelica.Mechanics.MultiBody.Frames.axesRotations(
  {1, 2, 3}, bodybox[1].body.angles_start, {0.0, 0.0, 0.0});

because in the latter case the left-hand-side has some indexing, while the former case doesn't.

I'm not sure whether this is because the new front-end creates a somewhat invalid data structure in this case, or rather that the back-end or simcode should be amended by adding the case of array element on the left-hand-side, maybe that's just a trivial fix.

@perost, @lochel, @adrpo, @wbraun can you agree on a solution for this case?

[Per Östlund]

Here's a small test model that gives the same issue when compiled using the old frontend:

model A
  parameter Modelica.Mechanics.MultiBody.Frames.Orientation R_start =
    Modelica.Mechanics.MultiBody.Frames.axesRotations(
      {1, 2, 3},
      angles_start,
      {0.0, 0.0, 0.0});

  parameter Real angles_start[3](each fixed = false);
initial equation
  angles_start = {0, 0, 0};
end A;

model test
  A a[1];
end test;

So it seems like this is something that should at least be fixed in the code generation.

[Francesco Casella]

As discussed at the devmeeting on Sep 17, @karim.abdelhak will have a look at this issue ASAP.

[Willi Braun]

Fixed with ​PR2678.

[Francesco Casella]

Unfortunately there are still issue with this model, as demonstrated e.g. by ​this report.

I made some further analysis on the optdaedump output. The flattened model contains the equation

revolute[2].frame_a.f = 
  -Modelica.Mechanics.MultiBody.Frames.resolve1(revolute[2].R_rel, revolute[2].frame_b.f)

which is then inlined and aliased by the back-end into

-bodybox[1].frame_b.f[1] = (-transpose(revolute[2].R_rel.T) * 
  {-bodyboxN.frame_a.f[1], -bodyboxN.frame_a.f[2], -bodyboxN.frame_a.f[3]})[1]
-bodybox[1].frame_b.f[2] = (-transpose(revolute[2].R_rel.T) *
  {-bodyboxN.frame_a.f[1], -bodyboxN.frame_a.f[2], -bodyboxN.frame_a.f[3]})[2]
-bodybox[1].frame_b.f[3] = (-transpose(revolute[2].R_rel.T) * 
  {-bodyboxN.frame_a.f[1], -bodyboxN.frame_a.f[2], -bodyboxN.frame_a.f[3]})[3] 

In the simulation results, bodyboxN.frame_a.f is nonzero due to the force applied on the beam tip, while bodybox[1].frame_b.f turns out to be zero, which is wrong, since the numerical value of revolute[2].R_rel.T turns out to be (correctly) close to the unit matrix.

I compared this with the result of optdaedump when using the old front-end, which is

-bodybox[1].frame_b.f[1] = (-revolute[2].R_rel.T[2,1]) * revolute[2].frame_b.f[2] - revolute[2].R_rel.T[3,1] * revolute[2].frame_b.f[3] - revolute[2].R_rel.T[1,1] * revolute[2].frame_b.f[1]
-bodybox[1].frame_b.f[2] = (-revolute[2].R_rel.T[2,2]) * revolute[2].frame_b.f[2] - revolute[2].R_rel.T[3,2] * revolute[2].frame_b.f[3] - revolute[2].R_rel.T[1,2] * revolute[2].frame_b.f[1]
-bodybox[1].frame_b.f[3] = (-revolute[2].R_rel.T[2,3]) * revolute[2].frame_b.f[2] - revolute[2].R_rel.T[3,3] * revolute[2].frame_b.f[3] - revolute[2].R_rel.T[1,3] * revolute[2].frame_b.f[1]

I don't know why the back-end expands the transpose operator after inlining the resolve1 function when using the OF, but it doesn't when using the NF. I guess this causes some issues with the generated code, which is not broken but gives wrong results. I flattened the model with both front-ends, and the two function printouts seem identical, so I don't really know where to look at.

@wbraun, @perost, any idea?

[Per Östlund]

Replying to casella:

I don't know why the back-end expands the transpose operator after inlining the resolve1 function when using the OF, but it doesn't when using the NF. I guess this causes some issues with the generated code, which is not broken but gives wrong results. I flattened the model with both front-ends, and the two function printouts seem identical, so I don't really know where to look at.

When using the OF it wasn't the backend that evaluated transpose, but the frontend. Or rather that the frontend expanded the argument to transpose and had a simplification rule for transpose of an array.

I tried to do the same in the NF, i.e. expand all function arguments, mostly to see what the impact would be. But the results were a bit mixed and on the whole didn't make much difference, so I disabled it. But I disabled it by just changing the default value of the debug flag I'd introduced, so -d=nfExpandFuncArgs can still be used to turn on expansion of all function arguments.

Instead I chose the more targeted approach of expanding the argument of only transpose in 5bbac02, which fixed the simulation of the FlexibleBeamModelica models as well as some of the MultiBody models in the MSL.

[Francesco Casella]

Replying to perost:

Instead I chose the more targeted approach of expanding the argument of only transpose in 5bbac02, which fixed the simulation of the FlexibleBeamModelica models as well as some of the MultiBody models in the MSL.

Excellent!