Error using variable coefficient in FENICS

Hello,

I am using FEATool to simulate a magnetostatic problem: a coil with an iron core. I implemented my own 2D axysimmetric equations (I also try the standard magnetostatic package), but I am having problems when using a variable dependent magnetic permeability (for the iron). I tried two approaches:
- Analytical formula describing magnetic permeability as function of B
- Tabulated data of magnetic permeability using an external function

In both cases the results are not congruent.
I read in another post that the FEATool solver might have trouble when using non-linear coefficients such as in this case, so I tried using the FENICS solver. The solver works well when using constant permeabilities but I have trouble when implementing a variable permeability. I tried both approaches:
- Analytical formula: it fails to adequately compute the dependent variable
- Tabulated data: seems to fail identifying the external function

I would appreciate any suggestions on how to approach this problem

Here are the two functions:
helmholtz2DFEA.m
mu_rel_iron.m

Thank you in advance

It seems to be a very non-linear problem, and I doubt in this case Fenics will manage better without some custom solution (you would also have to convert your external function and interpolation to Python).

So far the only approach I've found work is to manually linearize the "mu_rel" nonlinearity and iterate towards a solution while using this constant from the previous iteration, it is not very numerically efficient, but at least seems to converge:

helmholtz2DFEA.m
mu_rel_2.m

I'll update here if I can figure out a better or more efficient solution.

Thank you for the response!

I also manage to solve it using the tabulated function and reducing the non-linear relaxation parameter to 0.15-0.2, but it is not quite efficient.
I also want to try by transforming it into a time-dependent problem, and increase the current more smoothly. Maybe it could converge faster.

For the FENICS option, I also figure that I should convert my external function to python, however I was wondering where should I save the file to be accesible for FENICS? Or should I modify directly the converted python script?

Thank you again!

If your external function is small you should be able to just add it to the script. Otherwise put it so somewhere on the Python path, or add it as a sys import.

import sys
sys.path.append('/path/to/folder/with/python/files')
from my_python_file import my_function

After double checking, it seems there was an error in the previous linearization (evaluating the derivatives Ar and Az in mu_rel_2 to zero), and after that it doesn't converge anymore either. Moreover, a more traditional approach of slowly increasing the nonlinearity (and using the old solution as intial guess) like here does unfortunately not work either:

for nfac=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
  nfac

  fea.expr = { 'nfac', nfac;
               'Je', opt.Je;
               'mu_iron', 1;
               'mu0', 'pi*4e-7';
               'Br', '-Az';
               'Bz', 'A/r+Ar';
               'B', 'sqrt(Br^2+Bz^2)';
               'J', '0 0 Je';
               'mu_rel_2', '2890*exp(-((sqrt((A/r+nfac*Ar)^2+nfac*Az^2)-0.7)/0.66)^2)+1'
               'mu_rel', '1 mu_rel_2 1';
               'Bzz', '2*Ar' };

  fea.sol.u = solvestat( fea, 'nlrlx', 0.5, 'maxnit', 100, 'init', u0 );
  u0 = fea.sol.u;
end

helmholtz2DFEA2.m

Might worth a try with time stepping. I will update if I can figure out something better.

Thank you for the effort!

I forgot to attach in the first post the data for the external function:
BHcurve.mat

Also find attached my latest models:
helmholtz2DFEA.m
mu_rel_iron.m

I did some mesh refinement in the most relevant areas. In this way the computation time reduces even if using a small relaxation factor.
Also, I want to ask you if you have a better idea to adequately plot Bz. As you might have seen when r=0, Bz goes to Inf. I tried to handle this by approximating Bz using LHopital rule (what I called Bzz), and applying the following logical rule:

'Bz', '(r>0)*(A/r + Ar)+(r<=0)*2*Ar'

However, it seems that I am not avoiding the Inf, since it seems to be doing 0*Inf = Inf. Any suggestions?
Thank you again.

Yes, I checked the analytical problem as they should roughly have the same issues. Your new model looks good though. Regarding divisions of zero, or in your case nan (not-a-number from 0/0) two common approaches are:

1) Adding a very small number to the radius so that it will never be zero, but also not affected significantly for larger values " r -> r+sqrt(eps)"

2) For plotting you can use the setnan function (which returns the value if not nan, or zero if its nan) for example:

    postplot(fea,'surfexpr','setnan(B)')

If you are sure you have inf, you could make an equivalent function (in a plot nan show up as missing data, while inf would be some huge numbers):

function [ x ] = setinf( x, val )
if( isnumeric(x) )
  if( nargin<2 )
    val = 0;
  end
  ix = isinf(x);
  if( isscalar(val) )
    x(ix) = val;
  else
    x(ix) = val(ix);
  end
end

EDIT: The issue with expressions of type ("(r>0)*1/r + 2*(r<=0)") is that "(r>0)*1/r" will still be evaluated everywhere (switch expressions just evaluate to 0 or 1), so the result at r=0 will be 0*1/0 = nan.

Note that inputs to custom functions are not scalar (that would be far too slow), but arrays (usually in the quadrature points where the fem assembly occurs), which means that in your function "mu_rel_iron(r,A,Ar,Az)" the expression "if r < 0.1" is not scalar and I'm surprised it doesn't throw an error. Something like this should work better, although does not help with the convergence issues:

    B_magnitude = ((2.*Ar).^2+(-Az).^2).^(1/2);
    ix = r < 0.01;
    if( any(ix) )
      B_magnitude(ix) = ((A(ix)./r(ix) + Ar(ix)).^2+(-Az(ix)).^2).^(1/2);
    end