Laplace equation with integral constraint

Hello,
I would like to calculate the current distribution in a set of (superconducting) striplines (they are parallel and “infinite”).
Usually, the total current I_bias circulating in a (single) stripline is an input data/constraint for this problem (while the distribution in the stripline is far from being uniform). How would it be possible to implement such constraint ( iint from{Supercond} u (x,y)dS = I_bias ) with the basic installation of FreeFEM ?
or do I need to install additional software (IPOPT ?), and possibly recompile the code ?
Thank you for any advice.
Denis

I’m not sure I understand you correctly, but can’t you discretize your constraint using a varf in a vector c, and then solve this using a Lagrangian, e.g., matrix B = [[A,c],[c',0]];? Here is such an example where we impose a mean of 0. You can strip the PETSc stuff and use plain FreeFEM matrices.

Thank you very much for your suggestion. I have made quite a number a attempts to understand how to manipulate the stiffness matrix, the Lagrange function, and even the boundary conditions at the edge of the space domain. I may have to return to the basics and provide more details:

My objective is to calculate the current density in a system of superconducting wires. Input data are the geometry and the total currents circulating in the different wires. As in the general case the system does not have any symmetry (except translational invariance along direction z), I chose not to use cylindrical coordinates.

The script below simulates a single cylindrical infinite wire made with superconducting material and displays the following graphs: mesh, vector-potential component Az, and current density Jz in the wire.
For this particular example, the solution is well-known, Az ~ log® in the “air” (= outside the wire). Thus I can use the boundary condition Az=mu_0I/4pilog(x^2 +y^2) on the wall limiting the domain. Line 22 of the script can be modified to

  • change the space between rectangular section and circular sections of the box

  • change the script usng “PROBLEM/SOLVE” or using the stiffness matrix A : we obtain very similar results.
    Meissner-supra_wire-3.edp (4.7 KB)
    =====BEGINNING OF REPORT (Meissner-supra_wire-3.edp) ===============
    – FreeFem++ v 3.380001 (date Sat Feb 6 20:00:25 UTC 2016)
    Load: lg_fem lg_mesh lg_mesh3 eigenvalue
    1 : // Vector potential created by an infinite superconducting wire
    2 : // Solved using either the “problem” keyword, or manipulating the stiffness matrix.
    3 : // Boundary conditions at large distance being non-zero, a current must circulate
    4 : // Where Laplace(Az) is not 0, i.e. in the wire
    5 : // “boxoice” parameter controls whether we use a space with rectangular or cylindrical section.
    6 : ///////////////Define geometry///////////////////////////////
    7 : // linear dimensions (in µm)
    8 : // Penetration depth (squared)
    9 : real tgv=1E30, epsi=1E-6, xi,yi, Lond=0.2, l2= Lond^2; // Lambda_London =0.2 um
    10 : real mu0s2p = 2e-7, mu0=mu0s2p2pi, I1=1, A0=mu0s2pI1; // for current I1=1 Amp.
    11 : // Rectangular Section: Choice 1 = Results are marginally acceptable
    12 : real far=20,high=20;
    13 : border bottom(t=-far,far) { x=t; y=-high; label=1; }
    14 : border rside(t=-high,high) { x=far; y=t; label=1; }
    15 : border lside(t=high,-high) { x=-far; y=t; label=1; }
    16 : border top(t=far,-far) { x=t; y=high; label=1; }
    17 :
    18 : // Cylindrical section: Choice 2 = Results are just nonsense…
    19 : real RB=23; // RB is the radius of the box
    20 : border RoundBox(t=0,2
    pi) { x=RBcos(t); y=RBsin(t); label=2; }
    21 : Meissner-supra_wire_5bi4.edp
    22 : int boxoice= 2, prost=1; // boxoice: choice of outer box shape; prost: 0 for “problem”; 1 for “stiffness”
    23 :
    24 :
    25 : // Wire
    26 : real R=1; // R is the radius of the wire
    27 : border wire(t=0,2pi) { x=Rcos(t); y=Rsin(t); label=3; }
    28 : border bulk(t=0,2
    pi) { x=(R-2Lond)cos(t); y=(R-2Lond)sin(t); label=4; }
    29 : /////////////////////////////////////////////////////////////
    30 : //////////// Mesh and London equation ///////////////////////
    31 : int NBM=50, NBW=200;
    32 : mesh Th;
    33 : if ( boxoice==1 Meissner-supra_wire_5bi4.edpMeissner-supra_wire_5bi4.edp) Th=buildmesh(bottom(NBM)+rside(NBM)+top(NBM)+lside(NBM)+wire(NBW)+bulk(NBW));
    34 : if ( boxoice==2 ) Th=buildmesh(RoundBox(3
    NBM)+wire(NBW)+bulk(NBW));
    35 : // regions are numbered 0 , 1 and 2
    36 : int core=Th(0,0).region, super=Th(R-Lond,0).region, air=Th(R+Lond,0).region;
    37 : cout << “Superc. regions: " << core <<” and "<< super << ", and air: "<< air<< endl;
    38 :
    39 : fespace Vh(Th,P2); // With P1, the result would be very noisy…
    40 : Vh Az, v, u0, ubc;
    41 : int nbvertices=Th.nv, NDF=Vh.ndof;
    42 :
    43 : // The following lines give analytically exact BC for single cylindrical wire
    44 : if ( boxoice==1 ) {
    45 : u0 =(abs(abs(x)-far)<1E-2 || abs(abs(y)-high)<1E-2) ? A0/2
    log(x^2 +y^2):0; // for rectangular box
    46 : };
    47 : if ( boxoice==2 ) {
    48 : u0 =(abs(x^2+y^2-RB^2)<epsi) ? A0log(RB):0; // for disk box
    49 : };
    50 : ubc=(u0==0) ? 0: tgv
    u0; // This was to avoid 01E30, but it does not seem necessary …
    51 :
    52 : fespace W0(Th,P0);
    53 : W0 Sup=(region==super || region==core); // Size(Sup)=nb of TRIANGLES !!!
    54 : W0 Cor=(region==core);
    55 : cout << “Nb Vert.=”<< nbvertices<<", NDoF="<< NDF<< “, size(Sup)=” << Sup.n << endl;
    56 : func fL=-Sup/l2; // For London equation in the superconductor
    57 : ///////////////////////////////////////////////////
    58 : //////////// Solve problem ///////////////////////
    59 : problem poissonMagneto(Az,v) =
    60 : int2d(Th)( dx(Az)dx(v) + dy(Az)dy(v) )
    61 : + int2d(Th) ( -v
    Az
    fL ) // the source term: Mu0
    J=-Az/Lond^2 in the superconductor
    62 : + on(boxoice, Az=u0); // boundary condition far away: Az=mu0s2pI1log®
    63 :
    64 : varf va(Az,v) = // definition of the problem
    65 : int2d(Th)( dx(Az)dx(v) + dy(Az)dy(v))
    66 : + int2d(Th) ( -v
    Az
    fL ) // the source term: Mu0J=-Az/l2 in the superconductor
    67 : + on(boxoice,Az=u0);
    68 : matrix A=va(Vh,Vh);
    69 : real[int] F=ubc[];
    70 :
    71 : if (prost==0) {poissonMagneto;}
    72 : else {
    73 : set(A,solver=UMFPACK);
    74 : Az[]=A^-1
    F; // solve the linear system
    75 : }; // prost==1 seems to be slightly faster; otherwise, both methods give the same results.
    76 :
    77 : ////////////////////////////////////////////////////
    78 : //////////// Display results ///////////////////////
    79 :
    80 : plot(Th,ps=“wire.eps”,bw=1, wait=true, cmm=“Mesh”);
    81 :
    82 : Vh Bx = -dy(Az)(x,y), By = dx(Az)(x,y);
    83 : Vh Jz = (dx(By)-dy(Bx))/mu0;
    84 : string Title="Contour plot of Az, ";
    85 : if (prost==0) Title=Title+“Solve”; else Title=Title+“Stiffness”;
    86 : plot(Az, ps=“vector_potential_wire.eps”, value=true, wait=true, cmm=Title);
    87 : Title=“Contour plot of Jz, “;
    88 : if (prost==0) Title=Title+“Solve”; else Title=Title+“Stiffness”;
    89 : plot(Jz, ps=“Current_density_wire.eps”, value=true, wait=true, cmm=Title);
    90 : // Larger valuesMeissner-supra_wire_5bi4.edp of the current are indeed obtained for the surface of the wire,
    91 : cout << “Core current =” << int2d(Th)(CorJz) << endl; // 9% of the current in the core
    92 : // No negative values in the core, as expected
    93 : // Check total current: integrates the current density in the wire =I1
    94 : cout << “Total current =” << int2d(Th)(Jz) << " and Wire current =" << int2d(Th)(Sup
    Jz) << endl;
    95 : // The integrated current density is about 10% too small
    96 :
    97 : // Using Gauss theorem, calculate By(R,0)=mu0s2pI1/R, i.e. mu0s2pI1=RBy(R,0)
    98 : cout << "Apply Gauss theorem at r=2
    R: I1=”<< 2RBy(2R,0)/mu0s2p << endl;
    99 : // The value of the field is also too small (about -7%). Does not seem to depend on the mesh… ???
    100 : cout << "Apply Gauss theorem at r=20
    R: I1=”<< 20RBy(20*R,0)/mu0s2p << endl;
    101 : // Here also, By is too small (-9%).
    102 :
    103 : sizestack + 1024 =7872 ( 6848 )

    – mesh: Nb of Triangles = 42176, Nb of Vertices 21164
    Superc. regions: 0 and 1, and air: 2
    Nb Vert.=21164, NDoF=84503, size(Sup)=42176
    Core current =0.0295082
    Total current =83.7838 and Wire current =0.29396
    Apply Gauss theorem at r=2R: I1=0.302521
    Apply Gauss theorem at r=20
    R: I1=0.302628
    times: compile 0.018529s, execution 5.9965s, mpirank:0
    CodeAlloc : nb ptr 3284, size :368096 mpirank: 0
    Ok: Normal End
    ===========END OF REPORT==================
    This script gives unreliable results for circular section of the box limiting the space. When the space is defined with a rectangular section, the results are much better but not fully reliable:

  • Larger values of the current are indeed obtained for the surface of the wire,
  • No negative values in the core, as expected
  • The integrated current density is about 10% too small
  • the magnetic field “By” roughly follows a 1/r law.
  • The value of “By” is also too small (about -7%). Does not seem to depend on the mesh… ???

How is it possible to improve the results ?

The approaches presented above allow in principle to compute the current circulating in the wire. However, with a more complex system of conductors, the boundary conditions become difficult to evaluate beforehand. Is it possible to impose the value of the current I1 circulating in the wire (int2d(Th2, Jz) =I1) and calculate Az without giving boundary conditions on borders labeled “1” (or “2” in the case of the cylinbdrical box)?

I guess that it is necessary to manipulate the stiffness matrix A and the vector F associated to the linear form to implement the condition on the total currents circulating in the wires.
If I understand correctly, the equation A.u=F is not 100% correct because of edge effects: the points belonging to the border are missing neighbours to evaluate the 2nd derivative along the axis which is perpendicular to the border. This is presumably the fundamental reason why you must change the coefficients of the stiffness matrix for lines corresponding to points on the border, for a set of proper boundary conditions.

The second code (Meissner-supra_wire-5bi4.edp) attempts to implement the condition for the total current in a wire, using the algorithm of the Lagrange multiplier. This can be summarized as follows:
-find the solution of \Delta u = fLu given that \int u dS = mu0\London^2I, with fL=1/London^2 in the superconducting wire, and 0 elsewhere. For c = vector of triangle areas, it is equivalent to solve the matrix equations Au=F and c’u= II=mu0\Lambda^2*I, which result from the minimization of the Lagrage \int (U^2/2 -fL) dS - Lambda (\int u dS -I), where lambda is the Lagrange multiplier.
Therefore, I calculated the solution of AA * u= [[A, c];[c’,0]]U = [ZERO,II], where ZERO is the presumably the vector of NDoF elements, each of them equal to 0, and U = [u,Lambda] : Au+c
Lambda=ZERO; c’*u = II.

Because not giving boundary conditions is equivalent to giving Neumann boundary conditions, I still must include realistic BC, i.e. I must change the vector ZERO elements corresponding to the border. In a first try, I just want to give the Dirichlet boundary conditions corresponding to the expected analytical result, i.e. the same as in Meissner-supra_wire-3.edp.
The second code has 2 toggle parameters, boxoice and solvoice, the first one allows to choose between a rectangluar section box or a disk section box; the second parameter is used to choose either to solve the matrix exactly as in Meissner-supra_wire-3.edp, or to solve by minimization of the Lagrange function.
Meissner-supra_wire_5bi4.edp (6.7 KB)
========== BEGINNING OF REPORT (Meissner-supra_wire_5bi4.edp) =================
– FreeFem++ v 3.380001 (date Sat Feb 6 20:00:25 UTC 2016)
Load: lg_fem lg_mesh lg_mesh3 eigenvalue
1 : // Vector potential created by an infinite superconducting wire (circular cross-section)
2 : // Attempt to implement the total current condition with the “augmented matrix” approach
3 : // “boxoice” parameter controls whether we use a rectangular space or a disk.
4 :
5 : ///////////////Define geometry///////////////////////////////
6 : // linear dimensions (in µm)
7 : // Penetration depth (squared)
8 : real tgv=1E30, epsi=1E-10, xi,yi,Lond=0.2, l2= Lond^2; // Lambda_London =0.2 um
9 : real mu0s2p = 2e-7, mu0=mu0s2p2pi, I1=1, A0=mu0s2pI1; // for current I1=1 Amp.
10 :
11 : // Rectangular Box: Choice 1
12 : real far=20,high=20;
13 : border bottom(t=-far,far) { x=t; y=-high; label=1; }
14 : border rside(t=-high,high) { x=far; y=t; label=1; }
15 : border lside(t=high,-high) { x=-far; y=t; label=1; }
16 : border top(t=far,-far) { x=t; y=high; label=1; }
17 :
18 : // Disk box: Choice 2
19 : real RB=23; // RB is the radius of the box
20 : border RoundBox(t=0,2
pi) { x=RBcos(t); y=RBsin(t); label=2; }
21 :
22 : int boxoice= 2, solvoice=0; // boxoice: choice of outer box shape;
23 : // solvoice : choice of matrix. 0:modified stiffness matrix; 1:original stiffness matrix
24 :
25 : // Superconducting wire
26 : real R=1; // R is the radius of the superconducting wire
27 : border wire(t=0,2pi) { x=Rcos(t); y=Rsin(t); label=3; }
28 : border bulk(t=0,2
pi) { x=(R-2Lond)cos(t); y=(R-2Lond)sin(t); label=4; }
29 : /////////////////////////////////////////////////////////////
30 : //////////// Mesh and London equation ///////////////////////
31 : int NBM=50, NBW=250;
32 : mesh Th;
33 : if ( boxoice==1 ) Th=buildmesh(bottom(NBM)+rside(NBM)+top(NBM)+lside(NBM)+wire(NBW)+bulk(NBW));
34 : if ( boxoice==2 ) Th=buildmesh(RoundBox(3
NBM)+wire(NBW)+bulk(NBW));
35 : // regions are numbered 0 , 1 and 2
36 : int core=Th(0,0).region, super=Th(R-Lond,0).region, air=Th(R+Lond,0).region;
37 : cout << “Superc. regions: " << core <<” and "<< super << ", and air: "<< air<< endl;
38 :
39 : fespace Vh(Th,P2); // P1 gives too much noise.
40 : Vh Az, v, u0, ubc, Sup=(region==super || region==core); // Size(Sup)=Vh.ndof (nb of TRIANGLES with P0 and P1);
41 : int NDF=Vh.ndof;
42 : fespace W0(Th,P0);
43 : W0 Cor=(region==core);
44 : varf varea(unused,chiK) =int2d(Th)(chiK); // the space function constant / triangle
45 : real[int] c=varea(0,Vh);
46 : // I want to zero the elements of c not associated to a vertex (or triangle) in the superconductor
47 : // c=c.Sup[]; // Wrong: results do not look like results of Meissner-supra_wire-3.edp
48 : // I don’t understand why there is such a difference in integrating the current density only on the wire,
49 : // compared to integrating on the full space, i.e. on the wire + regions where the current density should be 0!!
50 : cout << “Size(Sup)=” << Sup.n << ", min© = " << c.min << ", max© = " << c.max << ", sum© = " << c.sum <<endl;
51 :
52 : func fL=-Sup/l2; // For London equation in the superconductor
53 :
54 : ///////////////////////////////////////////////////
55 : //////////// Solve problem ///////////////////////
56 :
57 : real lm;
58 : // Is “Not giving boundary conditions” equivalent to “Neumann” ?
59 : // It is OK for the case with cylindical symmetry, but not in the general case
60 : // In the rectangular space, do we need to give something more realistic ?
61 : // Or alternatively, how can we suppress the default Neumann condition ?
62 : // It should be enough to specify the source term with “I1”…
63 : // The following lines give analytically exact BC for single cylindrical wire
64 : if ( boxoice==1 ) {
65 : u0 =(abs(abs(x)-far)<1E-2 || abs(abs(y)-high)<1E-2) ? A0/2
log(x^2 +y^2):0; // for rectangular box
66 : ubc=(abs(abs(x)-far)<1E-2 || abs(abs(y)-high)<1E-2) ? tgv
A0/2log(x^2 +y^2):0; // for BC
67 : };
68 : if ( boxoice==2 ) {
69 : u0 =(abs(x^2+y^2-RB^2)<epsi) ? A0
log(RB):0; // for disk box
70 : ubc=(abs(x^2+y^2-RB^2)<epsi) ? tgvA0log(RB):0; // for BC
71 : };
72 : real[int] F = ubc[], FI=[F,-mu0l2I1], solu(NDF+1); // The source terms are all contained in FI +F (non0)
73 : cout << “NDoF=”<< NDF<<", size© = " << c.n << ", size(u0) = " << u0.n << endl;
74 :
75 : varf va(Az,vh) = // definition of the problem
76 : int2d(Th)( dx(Az)dx(vh) + dy(Az)dy(vh))
77 : + int2d(Th) ( -vh
Az
fL ) // the source term: Mu0J=-Az/l2 in the superconductor
78 : + on(boxoice,Az=u0) // applies the BC to the proper border
79 : ;
80 :
81 : matrix A=va(Vh,Vh), AA = [ [ A , c] , // As for Lagrange multiplier
82 : [ c’, 0] ] ; // “Total current” constraint: area’Az =Imu0s2p
83 :
84 : if (solvoice==0) {
85 : set(AA,solver=UMFPACK);
86 : solu=AA^-1
FI; // solve the linear system; BC are not correct at least for rectangular box
87 : // Instead of having iso-Az looking like circles, they gradually look like rectangles
88 : // when they approach the boundary of the domain, as if the whole boundary was set to the same value
89 : [Az[], lm] = solu;
90 : }
91 : else {
92 : set(A,solver=UMFPACK);
93 : Az[]=A^-1F; // solve the linear system: ONLY RECTANGULAR BOX gives results which are
94 : // reasonable = identical to what is found in Meissner-supra_wire-3.edp
95 : };
96 : ////////////////////////////////////////////////////
97 : //////////// Display results ///////////////////////
98 : // Retrieve the content of A and AA just for display
99 : int[int] I(1),J(1);
100 : real[int] C(1);
101 : [I,J,C]=A;
102 : cout <<“Size(A)=” << A.n << “, size(F)=” << F.n << “, Max(A)=”<< C.max<< “, Min(A)=”<< C.min<< endl;
103 : cout <<“Size(AA)=” << AA.n << “, Max(AA)=”<< C.max<< “, Min(AA)=”<< C.min<< endl;
104 : real MaxF=0, MinF=100, Fi;
105 : for (int i=0; i<F.n; i++) {
106 : if (F(i)>0) {
107 : Fi=F(i)/tgv;
108 : MaxF=MaxF>Fi? MaxF:Fi;
109 : MinF=MinF<Fi? MinF:Fi;
110 : };
111 : };
112 : cout <<"Max(F) = " << MaxF << “, meen(F)=”<< MinF<< “, tgv= “<< tgv<< endl;
113 : plot(Th,ps=“wire.eps”,bw=1, wait=true,cmm=“Mesh Box:”+NBM+” Wire:”+NBW);
114 :
115 : Vh Bx = -dy(Az)(x,y), By = dx(Az)(x,y);
116 : Vh Jz = (dx(By)-dy(Bx))/mu0;
117 : plot(Az, value=true, wait=true, cmm=“Contour plot of Az”);
118 :
119 : plot(Jz, value=true, wait=true, cmm=“Contour plot of Jz”);
120 : // Larger values of the current are indeed obtained for the surface of the wire,
121 : cout << “Core current =” << int2d(Th)(Cor
Jz) << endl; //
122 : // this is ~10% of the current in the surface, same ratio as in “Meissner-supra_wire-3.edp” :slight_smile:
123 : // No negative values in the core, as expected
124 :
125 : cout << "Lagrange factor : " << lm << “, size(Jz) = " << Jz.n << endl;
126 : // Check total current: integrates the current density in the wire =I1 ?
127 : cout << “Total current =” << int2d(Th)(Jz) << " and Wire current =” << int2d(Th)(SupJz) << endl;
128 : // => I do not get what I expect, the current in the wire not equal to I1
129 : // The integrated current density is by far too small. What is in the c matrix ???
130 :
131 : // Using Gauss theorem, calculate By(R,0)=mu0s2p
I1/R, i.e. mu0s2pI1=RBy(R,0)
132 : cout << "Apply Gauss theorem at r=2R: I1="<< 2RBy(2R,0)/mu0s2p << endl;
133 : cout << "Apply Gauss theorem at r=4R: I1="<< 4RBy(4R,0)/mu0s2p << endl;
134 : cout << "Apply Gauss theorem at r=20R: I1="<< 20RBy(20R,0)/mu0s2p << endl;
135 : // The above results show whether By(r,0) decreases as 1/r or not.
136 : // only in the central region of the domain for both rectangular and disk boxes
137 : // Is this problem related to the Neumann default boundary conditions ?
138 :
139 : sizestack + 1024 =11400 ( 10376 )

– mesh: Nb of Triangles = 53254, Nb of Vertices 26703
Superc. regions: 0 and 1, and air: 2
Size(Sup)=106659, min© = -1.14275e-16, max© = 0.300256, sum© = 1661.42
NDoF=106659, size© = 106659, size(u0) = 106659
– Block Matrix NxM = 2x2 nxm =106660x106660 nb none zero coef. 1429256
Size(A)=106659, size(F)=106659, Max(A)=1e+30, Min(A)=-2.48824
Size(AA)=106660, Max(AA)=1e+30, Min(AA)=-2.48824
Max(F) = 6.27099e-07, meen(F)=6.27099e-07, tgv= 1e+30
Core current =-0.0418992
Lagrange factor : 3.72496e-09, size(Jz) = 106659
Total current =88.6296 and Wire current =-0.412868
Apply Gauss theorem at r=2R: I1=-0.391594
Apply Gauss theorem at r=4
R: I1=-0.279846
Apply Gauss theorem at r=20*R: I1=3.2961
times: compile 0.018462s, execution 11.1552s, mpirank:0
CodeAlloc : nb ptr 3522, size :372400 mpirank: 0
Ok: Normal End
====================END OF REPORT===============

  • When the space is “rectangular” and solving for the “original matrix” A (boxoice= 1, solvoice=1), then the calculations are essentially identical to the cases (boxoice= 1, prost=0 or 1) in Meissner-supra_wire-3.edp. And the results are very comparable; it means that results are about 7…9% smaller than expected.
  • When the space is “rectangular” and solving for the “modified matrix” AA (boxoice= 1, solvoice=0), then the result become strange, as instead of having iso-Az looking like circles, they gradually look like rectangles when they approach the boundary of the domain, as if the whole boundary was set to the same value. I checked the values stored in F, and they look correct. Taking a look at the currents, the current density in the wire turned to NEGATIVE, which is against all expectations… and spikes of large positive currents show up at the surface of the wire. Total currents integrated over the different regions hardly depend on the value of the last element of FI (expected to fix the current in the wire).
  • When the space is cylindrical and solving for the “original matrix” A (boxoice= 2, solvoice=1), then the calculations are essentially identical to the cases (boxoice= 2, prost=0 or 1) in Meissner-supra_wire-3.edp. And the results are very comparable; it means that results are far off…
  • When the space is cylindrical and solving for the “modified matrix” AA (boxoice= 2, solvoice=0), then the results are as strange as in Meissner-supra_wire-3.edp; it means that results are far off, currents in the wire are NEGATIVE… etc. And the total current is not even comparable to the total current calculated with Meissner-supra_wire-3.edp, cylindrical cases.

I am convinced I missed something in the manipulation of the stiffness matrix, but I am really surprised to see that the shape of the space matters!

Thank you very much for any help.
Denis