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.

[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 https://github.com/modelica/ModelicaStandardLibrary/issues/1620 - and since the unitH in that code was incorrect I now more strongly believe that something else should be changed as indicated above.

[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 d0dbf4b1a8a3825dd5b84a5b01da3a8916fa8319 - 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?

[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.

[AHaumer]

To be honest: I'm not convinced. I'd leave eps without unit: parameter Real eps=mat.M/1000; and change the last two equations: hystR = -Js + unitT*(P1*P2-P3*P4) + mu0*Hstat - unitT*eps/2; hystF = Js + unitT*(P4*P2-P3*P1) + mu0*Hstat + unitT*eps/2; @christiankral what's your opinion?

[AHaumer]

Well even a better Solution: parameter SI.MagneticFluxDensity eps=unitT*mat.M/1000; and leave the two mentioned equations untouched. If that's not a killer, shouldn't we shift the milestone? @HansOlsson @henrikt-ma what's your opinion? I think we need more time to make all Hysteresis models consistent.

[HansOlsson]

see my last comment

Well, have you looked at #4241 ? As far I understand that presents how to make everything unit-correct in this model without any unitT - the only problem is that it has square root of units.

[henrikt-ma]

@HansOlsson @henrikt-ma what's your opinion?

I agree with @HansOlsson; this is just a short-term fix which we can comfortably apply, knowing that there's #4241 showing that there's something better we can do in the future (after adopting a future Modelica version where rational unit exponents are allowed).

[AHaumer]

@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.

[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 https://github.com/modelica/ModelicaStandardLibrary/pull/4111#issuecomment-1911686873.

[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.

[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. */