Machine overflow when defining boundary conditions

Question

Recently I have been trying to code Maxwell's equations over a closed surface and have been facing some trouble defining the boundary conditions for the magnetic field.

The equation for the normal of the magnetic field at the boundary is $B_1 = B_2$ where $B_1$ is on one side of the boundary, and $B_2$ is on the other side, and the distance between them is infinitesimally small.

Logically, we would then define the Dirichlet condition in NDSolve as such:

DirichletCondition[Bfx[x,y,z,t] == Bfx[x+10^-12,y,z,t], x== 0.5 || x== -0.5]

This code means that at the edges of my conductor, the $x$ component of the magnetic field a small distance apart should be equal.

However, this code would then give the error:

NDSolve::fembpw: The boundary condition {DirichletCondition[True,True]} cannot be parsed and will be ignored.

as well as the "arguments should be ordered consistently" error.

Some alternatives I have tried include equating the spatial derivatives of the magnetic field at the boundary to be zero, but since Maxwell's equations do not have such a high differential order, it will not be solved.

DirichletCondition[D[Bx[x,y,z,t],x]==0,x==0.5||x==-0.5]

Another method is Neumann values, but considering that I am solving for only three or six functions ($E$ field in $x$, $y$, $z$ and $B$ field in $x$, $y$, $z$), adding this extra Neumann value will cause the system to become over determined.

I would really appreciate it if someone who has experience in dealing with normal boundary conditions to give their advice.

Edit2 Now there are issues with zero pivot, which i found to be very hard to troubleshoot. Is there anything I should be looking out for? Here is the code

In[153]:= Needs["NDSolve`FEM`"];
w = 0.1;
h = 0.2;
th = 0.01;
\[Mu]0 = 4 \[Pi]*10^7;
\[Mu]c = 1.3*10^6;
m = {0, 0, 1};
\[Alpha] = 10^-6;
\[Beta] = 10^-6;
Pd[t_] := {Sin[t], Cos[t], 0};
\[Sigma] = 5.98*10^-7;
region = Cuboid[{-w/2, -w/2, h - 0.00001}, {w/2, w/2, h + th + 3}];
regionMesh = ToElementMesh[region];
Hm[x_, y_, z_, t_] := 
  1/(4 [Pi]) ((3 ( m . ({x, y, z} - Pd[t])) ({x, y, z} - Pd[t]))/
     Sqrt[(x - Sin[t])^2 + (y - Cos[t])^2 + z^2]^5 - m/
     Sqrt[(x - Sin[t])^2 + (y - Cos[t])^2 + z^2]^3);
mu[x_, y_, z_] := 
  If[-w/2 < x < w/2 || -w/2 < y < w/2 || 
    h < z < h + th, [Mu]c, [Mu]0] ;
Out[55]= (1.07907 z (x - 
   Sin[t]))/(z^2 + (y - Cos[t])^2 + (x - Sin[t])^2)^(5/2)
In[72]:= Monitor[sol = NDSolve[
   {
    Thread[
     Curl[{Efx[x, y, z, t], Efy[x, y, z, t], Efz[x, y, z, t]}, {x, y, 
        z}] == D[
        mu[x, y, 
          z](Hm[x, y, z, t] + {Hix[x, y, z, t], Hiy[x, y, z, t], 
            Hiz[x, y, z, t]}) , 
        t] + [Alpha] Laplacian[{Efx[x, y, z, t], Efy[x, y, z, t], 
          Efz[x, y, z, t]}, {x, y, z}]],
    Thread[
     Curl[
       mu[x, y, 
         z](Hm[x, y, z, t] + {Hix[x, y, z, t], Hiy[x, y, z, t], 
           Hiz[x, y, z, t]}), {x, y, 
        z}] == [Mu]0 [Sigma] {Efx[x, y, z, t], Efy[x, y, z, t], 
         Efz[x, y, z, t]} + [Alpha] Laplacian[{Hix[x, y, z, t], 
          Hiy[x, y, z, t], Hiz[x, y, z, t]}, {x, y, z}]],
    DirichletCondition[
     Hix[x, y, z, 
       t] == -Hm[x, y, z, t][[1]], (z == h || 
        z == h + th ) && (x == -w/2 || x == w/2)],
    DirichletCondition[
     Hiy[x, y, z, 
       t] == -Hm[x, y, z, t][[2]], (z == h || 
        z == h + th) && (y == -w/2 || y == w/2)],
    DirichletCondition[
     Hiz[x, y, z, t] == -Hm[x, y, z, t][[3]], (z == h || z == h + th)],
    DirichletCondition[
     Efx[x, y, z, t] == 
      0, ((z == h || z == h + th) && (y == -w/2 || y == w/2)) || (z ==
          h || z == h + th)],
    DirichletCondition[
     Efy[x, y, z, t] == 
      0, ((z == h || z == h + th) && (x == -w/2 || x == w/2)) || (z ==
          h || z == h + th)],
    DirichletCondition[
     Efz[x, y, z, t] == 
      0, ((z == h || z == h + th) && (y == -w/2 || 
          y == w/2)) || ((z == h || z == h + th) && (x == -w/2 || 
          x == w/2))],
    (*
    DirichletCondition[D[Hix[x,y,z,t],x]==0, h<z<h+th && (x==w/2 || 
    x==-w/2)],
    DirichletCondition[D[Hiy[x,y,z,t],y]==0, h<z<h+th && (y==w/2 || 
    y==-w/2)],
    DirichletCondition[D[Hiz[x,y,z,t],z]==0, z==h||z==h+th],*)
    Efx[x, y, z, 0] == 0,
    Efy[x, y, z, 0] == 0,
    Efz[x, y, z, 0] == 0,
    Hix[x, y, z, 0] == 0,
    Hiy[x, y, z, 0] == 0,
    Hiz[x, y, z, 0] == 0
    }, {Efx, Efy, Efz, Hix, Hiy, Hiz},
   {x, y, z} [Element] regionMesh, {t, 0, 1}, MaxStepSize -> 0.01,
   StepMonitor :> (p = t)
   ], p]

I tried to change up some conditions. So now my perpendicular magnetic field is zero at the boundaries and my tangential electric fields are zero at the boundaries as well. The conditions for the DirichletConditions may be a bit erroneous so I tried splitting it up into individual conditions but it didnt quite work either.

Now the main issue Im

Its a bit long but here goes:

DirichletCondition[
 Hix[x, y, z, 
   t] == \[Minus]Hm[x, y, z, t][[1]], (z == 
    h) && (x == \[Minus]w/2)],
DirichletCondition[
 Hix[x, y, z, 
   t] == \[Minus]Hm[x, y, z, t][[1]], (z == h) && (x == w/2)],
DirichletCondition[
 Hix[x, y, z, 
   t] == \[Minus]Hm[x, y, z, t][[1]], (z == 
    h + th) && (x == \[Minus]w/2)],
DirichletCondition[
 Hix[x, y, z, 
   t] == \[Minus]Hm[x, y, z, t][[1]], (z == h + th) && (x == w/2)],
DirichletCondition[
 Hiy[x, y, z, 
   t] == \[Minus]Hm[x, y, z, t][[2]], (z == 
    h) && (y == \[Minus]w/2)],
DirichletCondition[
 Hiy[x, y, z, 
   t] == \[Minus]Hm[x, y, z, t][[2]], (z == h) && (y == w/2)],
DirichletCondition[
 Hiy[x, y, z, 
   t] == \[Minus]Hm[x, y, z, t][[2]], (z == 
    h + th) && (y == \[Minus]w/2)],
DirichletCondition[
 Hiy[x, y, z, 
   t] == \[Minus]Hm[x, y, z, t][[2]], (z == h + th) && (y == w/2)],
DirichletCondition[
 Hiz[x, y, z, t] == \[Minus]Hm[x, y, z, t][[3]], (z == h)],
DirichletCondition[
 Hiz[x, y, z, t] == \[Minus]Hm[x, y, z, t][[3]], (z == h + th)],
DirichletCondition[
 Efx[x, y, z, t] == 0, (z == h) && (y == \[Minus]w/2)],
DirichletCondition[Efx[x, y, z, t] == 0, (z == h) && (y == w/2)],
DirichletCondition[
 Efx[x, y, z, t] == 0, (z == h + th) && (y == \[Minus]w/2)],
DirichletCondition[
 Efx[x, y, z, t] == 0, (z == h + th) && (y == w/2)],
DirichletCondition[
 Efy[x, y, z, t] == 0, (z == h) && (x == \[Minus]w/2)],
DirichletCondition[Efy[x, y, z, t] == 0, (z == h) && (x == w/2)],
DirichletCondition[
 Efy[x, y, z, t] == 0, (z == h + th) && (x == \[Minus]w/2)],
DirichletCondition[
 Efy[x, y, z, t] == 0, (z == h + th) && (x == w/2)],
DirichletCondition[
 Efz[x, y, z, t] == 0, (z == h) && (y == \[Minus]w/2)],
DirichletCondition[Efz[x, y, z, t] == 0, (z == h) && (y == w/2)],
DirichletCondition[
 Efz[x, y, z, t] == 0, (z == h + th) && (y == \[Minus]w/2)],
DirichletCondition[Efz[x, y, z, t] == 0, (z == h + th) && (y == w/2)],
```

Answer 1

When I tried to run your code, I got a series of messages indicating that many if not all of your DirichletConditions were being discarded. As was mentioned by Daniel Huber in the comments, it makes sense to use Neumann boundary conditions here. It's called NeumannValue instead of condition because you're supplying only the rhs of the condition. As explained in the documentation, the FEM only knows how to solve one equation - but it's a very general equation, which reduces to many common forms given appropriate choices for constants and various other expressions within. The Neumann boundary condition is obtained by doing an integration by parts on this equation, and the lhs take a predetermined form which you can find under NeumannValue/Details. What you need to supply is the flux across a surface along an outward going normal. In your case, there are no surface charges or surface currents, so both your E and H fields can be handled this way by specifying zero. This is the default behavior, so no boundary conditions have to be specified - NeumannValue[0,True] is automatically supplied. I.e. all fields normal to the boundary have the same value in and outside, so the derivatives are zero everywhere on the boundary.

There were a few other problems I found, or thought might be potential problems, so I've made a few changes. I'll detail them below.

`Needs["NDSolve`FEM`"];

I changed your code slightly, defining lower/upper limits for h, for convenience. I got rid of the small constant you subtracted at the lower limit, since I didn't see the point of it. You can restore it if you want. I also eliminated "th", since you seemed to be using that for you attempt to solve this with a set of Dirichlet conditions, which I eliminated in favor of Neumann conditions.

{w = .1, hllim=.2, hulim = (*3.21*).5, mu0 = 4*^7 Pi, muC = 1.3*^6, alpha = beta = 1*^-6, sigma = 5.98*^-7};

I use a Region here instead of an element mesh because it's easier to check if you're inside the region or not with \[Element]. I also removed it from NDSolve mainly to see if things ran any faster. You can restore it if you want; it can be added inside NDSolve, i.e. {x,y,z} \[Element] ToElementMesh[region] instead of what I have. I defined a function to make it easier to modify the size. I lost patience waiting for the code to complete, so I tried shrinking the size of the both the region and the time. Even with the current settings, it took ~3380 seconds to complete, with a region at least 10 times shorter and a time interval 10^5 shorter than yours.

constructRectangularRegion[llim_List, ulim_List] := 
  DiscretizeGraphics[Cuboid[llim, ulim]];
region = 
  constructRectangularRegion[{-w/2, -w/2, hllim}, {w/2, w/2, hulim}];

While trying to figure out what was the source of various problems, I ended up rewriting this. It should be equivalent to what you had.

Hm[x_, y_, z_, t_] :=
    With[{dCurr = {x, y, z} - {Sin[t], Cos[t], 0}, ez = {0, 0, 1}},
        mu[x, y, 
     z] ((3 ez . dCurr dCurr)/(dCurr . dCurr)^(5/2) - 
       ez/(dCurr . dCurr)^(3/2))/(4 Pi)];

I added an Attribute to mu so that differentiation treats it as a constant. I was concerned that you'd compute derivatives of it in your equations, and that NDSolve would start complaining about the number of equations being insufficient. I also used a test of a point being in/outside of the region instead of attempting to do this via inequalities. I think you had a mistake in your version: using Or instead of And.

SetAttributes[mu, Constant];
mu[pt__?NumericQ] /; {pt} \[Element] region := muC;
mu[__?NumericQ] := mu0;

For convenience, to save some typing.

eVec[x_, y_, z_, t_] := {Efx[x, y, z, t], Efy[x, y, z, t], 
   Efz[x, y, z, t]};
hVec[x_, y_, z_, t_] := {Hix[x, y, z, t], Hiy[x, y, z, t], 
   Hiz[x, y, z, t]};

As was the case for mu, I didn't want Hm to be treated as yet another set of differential equations to be solved, so I computed derivatives in the With, then inserted the lists inside the equations. The equations don't have Hm appear explicitly anywhere, just the computed values.

eqFaraday = With[{dHmt = D[Hm[x, y, z, t], t]},
    Thread[
    Curl[eVec[x, y, z, t], {x, y, z}] == 
     dHmt + D[hVec[x, y, z, t], t] + 
      alpha Laplacian[eVec[x, y, z, t], {x, y, z}]]];
eqAmpere = With[{curlHm = Curl[Hm[x, y, z, t], {x, y, z}]},
            Thread[curlHm + Curl[hVec[x, y, z, t], {x, y, z}] ==
 sigma eVec[x, y, z, t] + 
  alpha Laplacian[hVec[x, y, z, t], {x, y, z}]]];

ic = Join[Thread[eVec[x, y, z, 0] == 0], 
   Thread[hVec[x, y, z, 0] == Hm[x, y, z, 0]]];

As I mentioned above, I kept reducing the time range until I could get something to run. I'd meant to keep repeating this till I got something after 3 minutes, but got distracted and forgot to check. Almost an hour later, I got a result. The Method I added is recommended in the documentation. If you click on the red dots or blue circle when you get a message, it takes you to a page explaining the message. In this case, near the bottom of the page explaining the pivot error you saw, this recommendation was made. It seems to help.

sol = NDSolve[
  Flatten[{eqFaraday, eqAmpere, ic}], {Efx, Efy, Efz, Hix, Hiy, 
   Hiz}, {x, y, z} \[Element] region, {t, 0, 0.00001},
    Method -> {"IndexReduction" -> "StructuralMatrix"}]

I didn't get any error messages after this. Each InterpolatingFunction is roughly 9.7MB. Plotting was taking a long time, so I just dumped a few values with Table instead. Here's a screenshot.

A side note: you had some problems with the predicates you were using in DirichletCondition. The first couple returned true if satisfied your checks, and the x,y checks were superfluous. If you need to specify certain sub-regions on a boundary to operate on, you might want to consider using element markers. See the documentation forToElementMesh.

That's all I've got. I don't have a supercomputer handy, so I'm going to leave extending this to the full range and examining the results to you. Good luck.