# System of pde with Neumann boundary conditions

The Ginzburg-Landau equation for a system of squared superconductor with a slit in presence of a uniform magnetic field can be written as a sort of continuity equation with Neumann boundary conditions (Ref):

$$d_a. \frac{\partial \bf{u}}{\partial t} + \nabla. \Gamma = \bf{F}$$

$$\nabla \begin{pmatrix} u_1 \\ u_2 \end{pmatrix}.\bf{n}= 0$$

$$\nabla \times \begin{pmatrix} u_3 \\ u_4 \end{pmatrix} = \bf{B_a}$$

$$\begin{pmatrix} u_3 \\ u_4 \end{pmatrix}. \bf{n} = 0$$

The magnetic field is applied along the z-axis so $$\bf{B_a} = B_a \hat{k}, \bf{n} = \hat{k}$$, and $$B_a$$ is constant. The result is a set of four coupled pde :

dA = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, sigma, 0}, {0, 0, 0, sigma}};

u[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] := {u1[t, x, y, z], u2[t, x, y, z], u3[t, x, y, z],
u4[t, x, y, z]}

delgma[t_, x_, y_, z_] = {(-D[u1[t, x, y, z], {x, 2}] -
D[u1[t, x, y, z], {y, 2}])/(kappa^2), (-D[
u2[t, x, y, z], {x, 2}] -
D[u2[t, x, y, z], {y, 2}])/(kappa^2),
D[D[u4[t, x, y, z], x], y] - D[u3[t, x, y, z], {y, 2}] +
D[ba[x, y, z], y],
D[D[u3[t, x, y, z], x], y] - D[u4[t, x, y, z], {x, 2}] +
D[ba[x, y, z], x]};

f1[t_, x_, y_, z_] = -(D[u3[t, x, y, z], x] + D[u4[t, x, y, z], y])*
u1[t, x, y, z]/
kappa - (D[u1[t, x, y, z], x]*u3[t, x, y, z] +
D[u1[t, x, y, z], y]*u4[t, x, y, z])/
kappa - (u3[t, x, y, z]*u3[t, x, y, z] +
u4[t, x, y, z]*u4[t, x, y, z])*u1[t, x, y, z] + (1 - temp)*
u1[t, x, y,
z]*(1 - (u1[t, x, y, z]*u1[t, x, y, z] +
u2[t, x, y, z]*u2[t, x, y, z]));

f2[t_, x_, y_, z_] = -(D[u3[t, x, y, z], x] + D[u4[t, x, y, z], y])*
u2[t, x, y, z]/
kappa - (D[u2[t, x, y, z], x]*u3[t, x, y, z] +
D[u2[t, x, y, z], y]*u4[t, x, y, z])/
kappa - (u3[t, x, y, z]*u3[t, x, y, z] +
u4[t, x, y, z]*u4[t, x, y, z])*u2[t, x, y, z] + (1 - temp)*
u2[t, x, y,
z]*(1 - (u1[t, x, y, z]*u1[t, x, y, z] +
u2[t, x, y, z]*u2[t, x, y, z]));

f3[t_, x_, y_, z_] = (D[u2[t, x, y, z], x]*u1[t, x, y, z] -
D[u1[t, x, y, z], x]*u2[t, x, y, z])/
kappa - (u1[t, x, y, z]*u1[t, x, y, z] +
u2[t, x, y, z]*u2[t, x, y, z])*u3[t, x, y, z];

f4[t_, x_, y_, z_] = (D[u2[t, x, y, z], y]*u1[t, x, y, z] -
D[u1[t, x, y, z], y]*u2[t, x, y, z])/
kappa - (u1[t, x, y, z]*u1[t, x, y, z] +
u2[t, x, y, z]*u2[t, x, y, z])*u3[t, x, y, z];

f[t_, x_, y_, z_] = {f1[t, x, y, z], f2[t, x, y, z], f3[t, x, y, z],
f4[t, x, y, z]};

dudt[t_, x_, y_, z_] = {D[u1[t, x, y, z], t], D[u2[t, x, y, z], t],
D[u3[t, x, y, z], t], D[u4[t, x, y, z], t]};

eqns = dA.dudt[t, x, y, z] + delgma[t, x, y, z] - f[t, x, y, z];


The equations have to be solved over a squared region with a slit which I have tried to model as:

region3d =
ImplicitRegion[-1 <= x <= 1 && -1 <= y <=
1 && ! (0 < x <= 1 && -0.1 < y < 0.1) && 0.45 <= z <= 0.5, {x, y,
z} ];


I am trying to gave a shot at the solution with following parameters

ba[x_, y_, z_] = 1;
kappa = 4;
sigma = 2;
temp = 0.5;


I tried the following method and got the error:

Needs["NDSolveFEM"];

sol = NDSolveValue[{eqns[[1]] ==
NeumannValue[0, z == 0.45 || z == 0.5],
eqns[[2]] == NeumannValue[0, z == 0.45 || z == 0.5],
eqns[[3]] == NeumannValue[0, z == 0.45 || z == 0.5],
eqns[[4]] == NeumannValue[0, z == 0.45 || z == 0.5] ,
u3[0, x, y, z] == 0.5*x*ba[ x, y, z],
u4[0, x, y, z] == 0.5*x*ba[x, y, z]}, {u1, u2, u3, u4}, {t, 0,
10}, {x, y, z} \[Element] region3d,
Method -> {"FiniteElement",
"MeshOptions" -> {MaxCellMeasure -> 0.005}}];

NDSolveValue::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable. >>


Since I am solving such an equation for the first time, I am not particularly sure that I have used the Neumann boundary conditions correcly.

• Parameter dudt is not defined. Typo in boundary conditions: should be NeumannValue[0, z == 0.45 || z == 0.5]. Use {u1, u2, u3, u4} instead of u in NDSolve[]. Feb 7, 2020 at 14:20
• I've edited the question, I still get the same error after changing the BC.
– smj
Feb 7, 2020 at 18:47
• Now delete all [t_, x_, y_, z_]  from u, delgma, f, f1, f2,f3,f4,dudt,put f={f1,f2,f3,f4} and eqns = dA.dudt + delgma - f;. Add u1[0, x, y, z] == 0, u2[0, x, y, z] == 0, DirichletCondition[{u1[t, x, y, z] == 0, u2[t, x, y, z] == 0}, True], In NDSolve put eqns=={0,0,0,0}. Feb 7, 2020 at 22:14
• Using Method->"FiniteElement" is not ging to work for time dependent PDEs. Usw Method->{"MethodOfLines", .... } instead. You can find exampls on this site. Feb 8, 2020 at 5:35
• @ Alex Trounev I tried some of your suggestions but I still get the same error, moreover u1 and u2 do not satisfy any Dirichlet BC as shown in the equations I wrote above.
– smj
Feb 8, 2020 at 9:02

I propose to reproduce the 2D model from the article "Magnetic Flux Penetration in a Mesoscopic Superconductor with a Slit by Isaias G. de Oliveira". We use equations (1), (2) as follows

A = {u3[t, x, y], u4[t, x, y]};
psi = u1[t, x, y] + I u2[t, x, y];
sigma D[A, t] +
1/2/k I ((u1[t, x, y] -
I u2[t, x, y]) Grad[(u1[t, x, y] + I u2[t, x, y]), {x,
y}] - (u1[t, x, y] +
I u2[t, x, y]) Grad[(u1[t, x, y] - I u2[t, x, y]), {x,
y}]) + (u1[t, x, y]^2 + u2[t, x, y]^2) A -
Laplacian[A, {x, y}] // Simplify

D[psi, t] + I/k A.Grad[psi, {x, y}] - 1/k^2 Laplacian[psi, {x, y}] +
A.A psi -(1-T) psi (1 - (u1[t, x, y]^2 + u2[t, x, y]^2)) // Simplify

ComplexExpand[%]


Putting all the equations together we find

    eqns={-u1[t, x, y] + T*u1[t, x, y] + u1[t, x, y]^3 - T*u1[t, x, y]^3 + u1[t, x, y]*u2[t, x, y]^2 - T*u1[t, x, y]*u2[t, x, y]^2 + u1[t, x, y]*u3[t, x, y]^2 +
u1[t, x, y]*u4[t, x, y]^2 - (u4[t, x, y]*Derivative[0, 0, 1][u2][t, x, y])/k - Derivative[0, 0, 2][u1][t, x, y]/k^2 -
(u3[t, x, y]*Derivative[0, 1, 0][u2][t, x, y])/k - Derivative[0, 2, 0][u1][t, x, y]/k^2 + Derivative[1, 0, 0][u1][t, x, y],
-u2[t, x, y] + T*u2[t, x, y] + u1[t, x, y]^2*u2[t, x, y] - T*u1[t, x, y]^2*u2[t, x, y] + u2[t, x, y]^3 - T*u2[t, x, y]^3 + u2[t, x, y]*u3[t, x, y]^2 +
u2[t, x, y]*u4[t, x, y]^2 + (u4[t, x, y]*Derivative[0, 0, 1][u1][t, x, y])/k - Derivative[0, 0, 2][u2][t, x, y]/k^2 +
(u3[t, x, y]*Derivative[0, 1, 0][u1][t, x, y])/k - Derivative[0, 2, 0][u2][t, x, y]/k^2 + Derivative[1, 0, 0][u2][t, x, y],
u1[t, x, y]^2*u3[t, x, y] + u2[t, x, y]^2*u3[t, x, y] - Derivative[0, 0, 2][u3][t, x, y] + (u2[t, x, y]*Derivative[0, 1, 0][u1][t, x, y])/4 -
(u1[t, x, y]*Derivative[0, 1, 0][u2][t, x, y])/4 - Derivative[0, 2, 0][u3][t, x, y] + 2*Derivative[1, 0, 0][u3][t, x, y],
u1[t, x, y]^2*u4[t, x, y] + u2[t, x, y]^2*u4[t, x, y] + (u2[t, x, y]*Derivative[0, 0, 1][u1][t, x, y])/4 - (u1[t, x, y]*Derivative[0, 0, 1][u2][t, x, y])/4 -
Derivative[0, 0, 2][u4][t, x, y] - Derivative[0, 2, 0][u4][t, x, y] + 2*Derivative[1, 0, 0][u4][t, x, y]};


We use the data from the article and find a solution

L = 1; region2d =
ImplicitRegion[-L <= x <= L && -L <= y <=
L && ! (0 < x <= L && -0.1 < y < 0.1), {x, y}];
ba[x_, y_] := 1;
k = 4;
sigma = 2; T = .5; psi0 = 1/Sqrt[2 ];

ic = {u1[0, x, y] == psi0, u2[0, x, y] == psi0, u3[0, x, y] == 0,
u4[0, x, y] == 0};
bc = DirichletCondition[{u1[t, x, y] == psi0, u2[t, x, y] == psi0,
u3[t, x, y] == -0.5*y*ba[x, y] (1 - Exp[-5 t]),
u4[t, x, y] == 0.5*x*ba[x, y] (1 - Exp[-5 t])}, True];

t0 = 2; sol =
NDSolveValue[{eqns == {0, 0, 0, 0}, bc, ic}, {u1, u2, u3, u4}, {t, 0,
t0}, {x, y} \[Element] region2d];

Table[DensityPlot[sol[[i]][t0, x, y], {x, y} \[Element] region2d,
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotPoints -> 40, PlotRange -> All], {i, 4}]


I'm not sure that quantum vorticity can be obtained in this model. The author of the article missed the details. It is known that he used the COMSOL, but how the vortices arise there is not clear.

• It is a nice attempt, but the sign of the last two terms in D[psi, t] seems incorrect. Also I don't think u1 and u2 satisfy Dirichlet BC as you have mentioned. The equations seem to be soluble without these BC. Further to match the results like those shown in the paper, I think one has to take appropriate parameters such as L ~ 500 nm, t0 ~ 200 etc. Thanks.
– smj
Feb 9, 2020 at 18:59
• @smj You're right. I corrected the equation and calculated but did not receive the vortices. What does $L = 500 nm$ mean? In the article he took $L=5\lambda_0$. Feb 9, 2020 at 19:50
• Well $\lambda_0$ is the London penetration depth, which is typically around 100 nm. The author has not specified the exact values of $\lambda_0$ and $\sigma$ in the paper. Also I think the author has scaled $L$ and $t$ to dimensionless parameters.
– smj
Feb 10, 2020 at 16:31
• @smj After equations (1), (2) he wrote "In these equations, the distances are scaled according to the London penetration length λ and the time is scaled by ξ^2/D where D is a phenomenological diffusion coefﬁcient [12]." So in our model we can take L=5/2 and $\sigma=1$. Feb 10, 2020 at 17:05