Fix unit error in GenericHystTellinenEverett

[henrikt-ma]

Fixes unit error in binding equation:

parameter SI.MagneticFluxDensity eps=mat.M/1000;

[AHaumer]

Please add a description, and address the authors / library officers: @ThomasBoedrich and Johannes Ziske.

[christiankral]

Looks OK

[HansOlsson]

Note that it may be a different error instead (leading to unit errors); and based on the bad unit-casts in #4053 I think we need to investigate it more.

• One possibility would be that Modelica.Magnetic.FluxTubes.Material.HysteresisEverettParameter.BaseData should have SI.MagneticFluxDensity M instead of M(unit="1"); and then we should remove the other use of unitT in this model and in Modelica.Magnetic.FluxTubes.Utilities.everett and in Modelica.Magnetic.FluxTubes.Shapes.HysteresisAndMagnets.GenericHystPreisachEverett
• Another possibility would be that eps should relate to some other variable than mat.M; in the two preceding models it is Js and Br instead.
• It could, of course, be some more complicated error.

Even though unit-casting is sometimes needed, to me it smells too much like hiding the error without investigating it fully. But I don't have the detailed knowledge to investigate that (if I had copies of the references I could look into it).

However, to me it looks suspicious that after this change all uses of M will be multiplied by unitT.

For the units I could track down the discussion to #1620 - and since the unitH in that code was incorrect I now more strongly believe that something else should be changed as indicated above.

[HansOlsson]

Unit-checking is supposed to help us find and guide us towards correcting errors.

Here unit-checking has detected an error, and instead of understanding the error this PR just seems to silence the check-error. (My estimate is that 0.001 is just a multiplicative small factor in this example - similar to the other similar models; I don't see that it means multiplying the unit-less number M with 0.001 Tesla as some universal small magnetic field - (about 3 times earth's magnetic field)).

That is not how we are supposed to use unit-checking.
It could instead be that it's a mistake in the formula, that M has the wrong unit, or something else.

[henrikt-ma]

That is not how we are supposed to use unit-checking. It could instead be that it's a mistake in the formula, that M has the wrong unit, or something else.

You are right. The intention of the PR is to fix the error, but since I am not an expert in this domain I am not surprised that my initially suggested "fix" isn't deemed the right one. If a domain expert steps forward and suggests the proper way of fixing the unit error, I'll happily go with that fix instead of mine.

[henrikt-ma]

@HansOlsson, as far as I can tell, this change is in the same spirit as the other PRs that fix unit errors, which the previous MAP-Lib meeting decided to go ahead with. Would you please unblock?

[HansOlsson]

@HansOlsson, as far as I can tell, this change is in the same spirit as the other PRs that fix unit errors, which the previous MAP-Lib meeting decided to go ahead with. Would you please unblock?

As previously stated it's not:

This model involves a material with

   parameter Real M(final unit="1")=0.95
    "Related to saturation value of magnetization";
  parameter Real r(final unit="1")=0.55
    "Proportion of the straight region in the vicinity of Hc";

and also an equation with unit-casting inside it: hystR = -Js + unitT*(P1*P2-P3*P4) + mu0*Hstat - eps/2;, and parameter Real P1 = (mat.M*mat.r*(....

So we already have one unit-cast, and several variables with unit="1"; combined this indicates some other underlying unit-error(s), and this PR is just trying to patch it over without investigating the actual error.

I don't know what the actual unit-error is, in contrast to #4053 and #4158 (for related models) where it was fairly straightforward to find the error; but I doubt that all of M, P1, P2 should really have unit="1".

Trying M.unit="T" gives:

parameter SI.MagneticFluxDensity eps=mat.M/1000;
hystR = -Js + (P1*P2-P3*P4)/unitT + mu0*Hstat - eps/2;

and P1.unit="T", P2.unit="T"; which is consistent with only one unit-cast - but it's still not perfect.

It could also be that there's some other error.

[henrikt-ma]

So we already have one unit-cast, and several variables with unit="1"; combined this indicates some other underlying unit-error(s), and this PR is just trying to patch it over without investigating the actual error.

Why don't we do both:

• This PR to make the model pass fundamental unit checks.
• Report an issue about overall dubious use of units.

[HansOlsson]

So we already have one unit-cast, and several variables with unit="1"; combined this indicates some other underlying unit-error(s), and this PR is just trying to patch it over without investigating the actual error.

Why don't we do both:

• This PR to make the model pass fundamental unit checks.
• Report an issue about overall dubious use of units.

Because it doesn't fix anything, and it might be the wrong change.

It's like adding pointer-casts in C code instead of looking up what the correct arguments should be and using them; basically silencing the compiler instead of improving the model - which makes it more likely that actual issue go uncorrected.

It's just a distraction from finding the actual issue.

Someone needs to look up what the units should be in this model, and preferably add a reference that explains them.

[HansOlsson]

See #4241 for an attempt to actually clean up these models.
It was split into different commits so that it might be possible to cherry-pick the ones that work in Modelica Language 3.6 (not guaranteed to work, it might also create a mess).

Note that there might be errors in the equations (instead of/in addition to) these unit-errors. I simply don't know, and I couldn't find the information in the references.

[henrikt-ma]

See #4241 for an attempt to actually clean up these models. It was split into different commits so that it might be possible to cherry-pick the ones that work in Modelica Language 3.6 (not guaranteed to work, it might also create a mess).

With #4241 representing the bright future when MSL uses some Modelica specification which hasn't been released yet, shouldn't we avoid the potential mess by just hacking the equations in the simplest way (represented by the this PR) to make them unit consistent in the meantime?

[HansOlsson]

See #4241 for an attempt to actually clean up these models. It was split into different commits so that it might be possible to cherry-pick the ones that work in Modelica Language 3.6 (not guaranteed to work, it might also create a mess).

With #4241 representing the bright future when MSL uses some Modelica specification which hasn't been released yet, shouldn't we avoid the potential mess by just hacking the equations in the simplest way (represented by the this PR) to make them unit consistent in the meantime?

One conclusion in #4241 is that the line that is equation that this PR tries to unit-correct is just plain wrong (but it is just an epsilon so it doesn't matter that much), and thus the proposal is to replace it with another equation in d0dbf4b - and that specific change doesn't require any new features.

That PR represent what unit-checking is supposed to do: help us find errors in equations.
Not hiding unit-errors by sprinkling unit* in equations to pass some test.

[henrikt-ma]

One conclusion in #4241 is that the line that is equation that this PR tries to unit-correct is just plain wrong (but it is just an epsilon so it doesn't matter that much), and thus the proposal is to replace it with another equation in d0dbf4b - and that specific change doesn't require any new features.

Well, that commit by itself doesn't look like a mess to me. Could you please turn that into a separate PR that would replace this one?

[casella]

@ThomasBoedrich resigned as library officer two weeks ago, so the ultimate decision about this PR rests with the new library officers @AHaumer and @christiankral.

@christiankral already approved this. @AHaumer what do you think?

[AHaumer]

ok for me

[HansOlsson]

One conclusion in #4241 is that the line that is equation that this PR tries to unit-correct is just plain wrong (but it is just an epsilon so it doesn't matter that much), and thus the proposal is to replace it with another equation in d0dbf4b - and that specific change doesn't require any new features.

Well, that commit by itself doesn't look like a mess to me. Could you please turn that into a separate PR that would replace this one?

To be clear the line is: parameter SI.MagneticFluxDensity eps=Js/1000;

As far as I understand the slight change in eps doesn't matter in practice.

Yes, I could do that as soon as someone could review #4241 and see that it long-term is the correct solution (avoiding unit-casting completely in the models, but instead having square root of T as unit which will have to wait for Modelica Language 3.7).

[henrikt-ma]

Requesting new reviews, since the recent change from @HansOlsson changes the PR entirely.

[AHaumer]

To be honest: I need some time to dive into that topic.

[HansOlsson]

Ok now.

[AHaumer]

For me it seems to be a good short-term fix until taking care of #4241

[henrikt-ma]

For me it seems to be a good short-term fix until taking care of #4241

Good. @AHaumer then you should be able to merge this, and ask @Harisankar-Allimangalath to also include the fix for 4.1.0.

Edit: See #4111 (comment).

[HansOlsson]

@HansOlsson yes I've looked at #4241, we will care about that later (for 4.2.0). Yes here we could use a short-term fix.

But this short-term fix is just saying: "The equation (and units all over the model) are wrong here as far as we can see, but instead of fixing that we instead make the unit-checker shut up by multiplying with unitT."

[henrikt-ma]

@HansOlsson yes I've looked at #4241, we will care about that later (for 4.2.0). Yes here we could use a short-term fix.

But this short-term fix is just saying: "The equation (and units all over the model) are wrong here as far as we can see, but instead of fixing that we instead make the unit-checker shut up by multiplying with unitT."

Sorry @HansOlsson, I missed that @AHaumer pushed 6df5d99d. (I wasn't aware that others could push to my fork like this.) I take back my suggestion to merge.

[AHaumer]

I'm against the solution changing the value of eps, this influences the results.
It seems that we don't find a solution to which anyone agrees, therefore I'm shifting the milestone for a solution together with #4241

[AHaumer]

@HansOlsson I don't understand why equation and units all over the model should be wrong?
Ok we have to check all parameters and equations of all hysteresis models.
Unfortuantely a lot of parameters are defined as Real: HysteresisEverettParameter.BaseData.{M, r, q, p1, p2, K, sigma} without proper explanation.
parameter SI.MagneticFluxDensity eps=mat.M/1000; is used in 2 equations:
hystR = -Js + unitT*(P1*P2-P3*P4) + mu0*Hstat - eps/2;
hystF = Js + unitT*(P4*P2-P3*P1) + mu0*Hstat + eps/2;
It is wrong (and not backwards compatible) to replace mat.M/1000 by Js/1000.
What's the problem with the equations?
hystR and hystF are MagneticFluxDensity, as well as Js.
Not very neat unitT*(P1*P2-P3*P4) is MagneticFluxDensity.
mu0*Hstat is MagneticFluxDensity, too.
Ok I admit 'parameter Real mu0(final unit="N/A2") = mu_0;' is not 100% clean.
Here Modelica.Units.SI.Permeability should be used: N/A2 = H/m = V.s/A.m
Last term: eps is MagneticFluxDensity.
@HansOlsson and @henrikt-ma could you explain the problems and what's your way to solve them?

[henrikt-ma]

@HansOlsson and @henrikt-ma could you explain the problems and what's your way to solve them?

Now that I see that you are speaking in the role of the library officer, I see things differently. As the current state of the PR is the same as the original quick-fix by me, I was afraid of pushing for this to be merged after there had been a discussion resulting in a different solution. Now, if you as the library officer see a problem with that solution, I have no reason to object to going back to the original short-term quick-fix.

[HansOlsson]

@HansOlsson I don't understand why equation and units all over the model should be wrong? Ok we have to check all parameters and equations of all hysteresis models. Unfortuantely a lot of parameters are defined as Real: HysteresisEverettParameter.BaseData.{M, r, q, p1, p2, K, sigma} without proper explanation. parameter SI.MagneticFluxDensity eps=mat.M/1000; is used in 2 equations: hystR = -Js + unitT*(P1*P2-P3*P4) + mu0*Hstat - eps/2; hystF = Js + unitT*(P4*P2-P3*P1) + mu0*Hstat + eps/2;

Let me try to explain:

I agree that a lot of variables are declared as Real without any proper explanation, and we don't really have a good explanation for some of the units for the variables that have unit="1" either (like M). That is a minor warning sign.

Additionally we have several multiplications with unitT. This is a big warning sign.

Multiplying by such an unit-factor makes sense in s-parametrizations (where we don't really compare the different values), if the value are part of an external interface (where units are missing), and in some other odd cases - but I cannot see that it makes sense here. And #4053 show that related models got it wrong.

To me this indicates that the multiplication by unitT is hiding an underlying unit-error, and we need to understand what it really is (and if possible remove it) - and to me the conclusion in #4241 is that:
P1.unit="T(1/2)" (meaning square root of Tesla), mat.r.unit="1", and mat.M.unit="T(1/2)" (instead of "1").
That allow us to remove the multiplications by unitT, in three different classes; which to me is a strong indication that there was a general issue with the units here. However, that unit-syntax isn't in Modelica 3.6.

The only thing that doesn't make sense with those units is the eps-formula, but I believe that the previous model was just plain wrong - (in an un-important way) and eps was just intended to be a small value, and the exact value really doesn't seem to matter (for the example models).

It may be that the conclusion was wrong, and that the units should be something else, or there is some other issue else hiding, - but my point remain that we need to understand the formulas, and what the units should be, and not just multiply by unit* to satisfy some unit-checking - instead we should use unit-checking to discover these issues; and propose possible solutions.

[HansOlsson]

This just hides the underlying unit-errors as explained in the longer explanation - at least as far as I can tell.

[henrikt-ma]

This just hides the underlying unit-errors as explained in the longer explanation - at least as far as I can tell.

I understand this point, but I think of it as this PR changing a bad model with unit inconsistencies into a bad model with consistent units. Regardless of this change, we have #4241 as a reminder that the model needs improvement, but in the meantime we could at least have a model that shouldn't fail to build for users that opt in for treating unit inconsistencies as errors (either according to a future specification, or according to tool-specific rules if there is no specification to align with).

[HansOlsson]

This just hides the underlying unit-errors as explained in the longer explanation - at least as far as I can tell.

I understand this point, but I think of it as this PR changing a bad model with unit inconsistencies into a bad model with consistent units.

But the "consistent units" are, as far as I can tell, based on incorrect units and I would only describe it as "passing unit checks" - not as having "consistent units". This mean that there's a risk that those bad units will spread more; and after correcting the underlying unit errors we will have to update this equation once more.

[henrikt-ma]

Maybe we could look into a workaround for #4241 while waiting for rational exponents in the unit expression syntax. We could try the following pattern:

  parameter Real p = 42; /* Will get inferred unit T(1/2) */
protected
  final parameter Real p2(unit = "T") = p^2; /* Workaround for directly setting unit = "T(1/2) on p. */

[AHaumer]

@henrikt-ma smart idea. I tried something similar with Dymola 2024x but without success:

  constant Modelica.Units.SI.MagneticFluxDensity B=sqrtB^2;
  constant Real sqrtB=1;

I suppose we have to be patient:

• wait for rational exponents syntax
• go through all hysteresis models, write together all constants, parameters and equations (and function calls).
• correct the units of all parameters

[HansOlsson]

@henrikt-ma smart idea. I tried something similar with Dymola 2024x but without success:

  constant Modelica.Units.SI.MagneticFluxDensity B=sqrtB^2;
  constant Real sqrtB=1;

I suppose we have to be patient:

• wait for rational exponents syntax
• go through all hysteresis models, write together all constants, parameters and equations (and function calls).
• correct the units of all parameters

You can get it to work, sort of, in Dymola:

parameter Real p=42;
protected 
 final parameter Real p2(unit="T", fixed=false);
initial equation 
 sqrt(p2)=p;

when activating: Advanced.Modelica.CheckUnitsProposedSimple = true;

We will see if we can the other variant to work for the next release. However, for constants it may be slightly messier.

[henrikt-ma]

Yes, the idea doesn't work for constants, at least not if they are handled according to current plans for standardization.

Actually, I shouldn't even have proposed this as a workaround with parameters, as I consider it quite problematic to allow expressions with empty unit to be implicitly associated with a unit determined by unit inference. In my opinion the way to do it would require changing the parameterization so that the user specifies a non-negative value in tesla, while the original parameter is made finally equal to the square root. Unfortunately, this would be a backwards incompatible change, as code modifying the original parameter would become illegal.

[AHaumer]

@HansOlsson how could I check if this weird unit-cast is succesfull?
Slightly different to your test:

  parameter Real sqrtB=1;
protected 
  final parameter Modelica.Units.SI.MagneticFluxDensity B(fixed=false);
initial equation 
  B=sqrtB^2;

@henrikt-ma yes you're right: This is a dirty hack, it's much better to declare the units explicitly.
In the meantime I'll check the Hysteresis models and parameters.
I suppose this leads to understand the modeling of hysteresis. Unfortunately the documentation is very useful.

see later discussions

[AHaumer]

There are three type of hysteresis models:

• TellinenSoft and TellinenHard contain the parameter definitions in the model
• TellinenEverett and PreisachEverett (based on the Everett-function) use the parameter record from HysteresesEverettParameter
• TellinenTable use the parameter record HysteresisTableData

Regarding the tables:
The define the rising and the falling branch separately.
Each branch is defined by a Real[:,2]-table: first column is magnetic field strength H, second column is flux density B.
Is there any way to define the columns of a table with different units?

[henrikt-ma]

Each branch is defined by a Real[:,2]-table: first column is magnetic field strength H, second column is flux density B.
Is there any way to define the columns of a table with different units?

I think this is a topic way ahead of the current discussions around unit checking. Ideally, would have liked this to work:

parameter Real[:, 2] table(each unit = {"A/m", "T"}) = …;

However I'm afraid the definition of each does not allow this interpretation of a multi-dimensional array as a nested array.

A workaround would be to define the two columns separately:

parameter Real[:] fieldStrength(each unit = "A/m") = {1, 2, 3};
parameter Real[:] fluxDensity(each unit = "T") = {4, 5, 6};
final Real[:, 2] table = [fieldStrength, fluxDensity];

Of course, this makes it very inconvenient to edit the H-B data as sequence of pairs.

[HansOlsson]

@HansOlsson how could I check if this weird unit-cast is succesfull?

Unfortunately it seems that support was added after Dymola 2024x was released.
But generally you translate the model and look in the plot-browser (or plot the variable); you should see weird units for these variables.

Slightly different to your test:

  parameter Real sqrtB=1;
protected 
  final parameter Modelica.Units.SI.MagneticFluxDensity B(fixed=false);
initial equation 
  B=sqrtB^2;

That required additional logic - since previously Dymola didn't deduce the unit of a in a^2=b; - it only propagated the unit of a forward to b. (It will be improved.)

[AHaumer]

@HansOlsson and @henrikt-ma I have some questions:

• Are units checked within functions?
• Is it checked whether the function's return value fits into the equation where the function is called?
• Should we shift the discussion to Hysteresis37 #4241 ?

I've already started to check the units in all hysteresis models. As far as I see:

• The partial models GenericHysteresis and GenericHysteresisTellinen are ok (but they influence the models that extend from them)
• GenericHystTellinenSoft is ok
• GenericHystTellinenHard seems to be buggy:

parameter Real M = 10/Hc "Slope of tanh()-function";
final parameter SI.MagneticFieldStrength H0=0.5*log((1 + mu0*Hc/Br)/(1 - mu0*Hc/Br)) + M*Hc;;

• HysteresisTableData.BaseData and GenericHystTellinenTable seem to be ok.
• HysteresisEverettParameter.BaseData, function everett, function initPreisach, GenericHystTellinenEverett, GenericHystPreisachEverett are hard to read - it will take some time to analyze them.
  I'll also try to get [YUY89] | Yamaguchi, T.; Ueda, F.; Yamamoto, E.: Simulation of Hysteresis Characteristics of Core Materials Using the Everett Function, IEEE Translation Journal on Magnetics in Japan, vol.4, no.6, pp.353-359, 1989.
  At the end, I'd like to get ALL hysteresis models clean.

[henrikt-ma]

@HansOlsson and @henrikt-ma I have some questions:

• Are units checked within functions?
• Is it checked whether the function's return value fits into the equation where the function is called?

There can only be tool-specific answers at the moment, plus anecdotal mentions of various ideas for standardization being discussed – both within and outside the language group.

At Wolfram, we are convinced that functions shouldn't escape unit checking, and the basic principles for how it could be defined see to follow naturally from the general principles of our unit checking machinery. We are working on a proposal for standardized unit checking based on the general principles of our machinery, but as we are aware of the poor status of unit consistency inside existing functions we plan to not propose that function bodies are considered in a first version of standardized unit checking. I hope, however, that there will still be tool-specific ways to opt in for also considering the bodies of functions.

In the proposal we are working on, declared units on function inputs and outputs will be considered.

• Should we shift the discussion to Hysteresis37 #4241 ?

Yes, the longer discussion about how to do it right would fit better in #4241; here we should just decide about what to do short-term.

I've already started to check the units in all hysteresis models. As far as I see:

It is great to hear that you have initiated this work! Hopefully it will not be that complicated to also complete it – unless we hit the road block of not being able to make needed backwards-incompatible changes, such as replacing a "T(1/2)" parameter by a "T" parameter…

[HansOlsson]

@HansOlsson and @henrikt-ma I have some questions:

• Are units checked within functions?

Dymola has the possibility, but as default doesn't.

• Is it checked whether the function's return value fits into the equation where the function is called?

Dymola normally checks that. However, there were built-in functions that slipped by earlier - see below.

Regarding functions (and classes in general) it is important to note that each call of the function is checked somewhat separately (so the same function can add "kg", "A" or "m" as long as they are different calls).

• Should we shift the discussion to Hysteresis37 #4241 ?

I've already started to check the units in all hysteresis models. As far as I see:

• The partial models GenericHysteresis and GenericHysteresisTellinen are ok (but they influence the models that extend from them)
• GenericHystTellinenSoft is ok
• GenericHystTellinenHard seems to be buggy:

parameter Real M = 10/Hc "Slope of tanh()-function";
final parameter SI.MagneticFieldStrength H0=0.5*log((1 + mu0*Hc/Br)/(1 - mu0*Hc/Br)) + M*Hc;;

That was already corrected in #4053 replacing the units in those declarations.
The declaration above (from before the correction) are a complete mess regarding units - and the correction included removing unitH-casting among the equations.

At the end, I'd like to get ALL hysteresis models clean.

I agree that would be good.