# Implementation of MQS equations and results

## 1. Implementation

Here are the important points of the code that solves the MQS system of equations in our test cases. The file is called "mqs-form.cpp" in the master branch.

We first define our function spaces Ah and Vh, as well as A and V on their respective space :

``````  auto Ah = Pchv<1>( mesh );
auto Vh = Pch<1>( mesh, markedelements(mesh, range) );

auto A = Ah->elementPtr();
auto V = Vh->elementPtr();``````

We then initialize A and V by A0 and V0 given :

``````  auto A0 = expr<3, 1>(soption(_name="A0"));
auto V0 = expr(soption(_name="V0"));
(*A) = project(_space = Ah, _expr = A0);
(*V) = project(_space = Vh, _expr = V0);``````

We then define our "blocks" where we will store our equations. The block structure consists of transforming our system of equations into the form $MU=F$ where U is the desired vector.

``````  BlocksBaseGraphCSR myblockGraph(2,2);
myblockGraph(0,0) = stencil(_test=Ah,_trial=Ah, _diag_is_nonzero=false, _close=false)->graph();
myblockGraph(0,1) = stencil(_test=Ah,_trial=Vh, _diag_is_nonzero=false, _close=false)->graph();
myblockGraph(1,0) = stencil(_test=Vh,_trial=Ah, _diag_is_nonzero=false, _close=false)->graph();
myblockGraph(1,1) = stencil(_test=Vh,_trial=Vh, _diag_is_nonzero=false, _close=false)->graph();
auto M = backend()->newBlockMatrix(_block=myblockGraph);

BlocksBaseVector<double> myblockVec(2);
myblockVec(0,0) = backend()->newVector( Ah );
myblockVec(1,0) = backend()->newVector( Vh );
auto F = backend()->newBlockVector(_block=myblockVec, _copy_values=false);

BlocksBaseVector<double> myblockVecSol(2);
myblockVecSol(0,0) = A;
myblockVecSol(1,0) = V;
auto U = backend()->newBlockVector(_block=myblockVecSol, _copy_values=false);``````

We then implement our equations as well as the boundary conditions :

``````  for (t = dt; t < tmax; t += dt)
{
auto M00 = form2( _trial=Ah, _test=Ah ,_matrix=M, _rowstart=0, _colstart=0 );
for( auto const& pairMat : M_modelProps->materials() )
{
auto name = pairMat.first;
auto material = pairMat.second;
auto mur = material.getScalar("mu_mag");

M00 += integrate( _range=markedelements(mesh, material.meshMarkers()),
}

auto M01 = form2( _trial=Vh, _test=Ah ,_matrix=M, _rowstart=0, _colstart=1 );
auto F0 = form1( _test=Ah, _vector=F, _rowstart=0 );
auto M11 = form2( _trial=Vh, _test=Vh ,_matrix=M, _rowstart=1, _colstart=1 );
auto M10 = form2( _trial=Ah, _test=Vh ,_matrix=M, _rowstart=1, _colstart=0 );
auto F1 = form1( _test=Vh ,_vector=F, _rowstart=1 );

for( auto const& pairMat : M_materials )
{
auto name = pairMat.first;
auto material = pairMat.second;
auto sigma = material.getScalar("sigma");

M00 += integrate(_range=markedelements(mesh, material.meshMarkers()),
_expr = mu0 * sigma * inner(id(A) , idt(A) ));

M01  += integrate(_range=markedelements(mesh, material.meshMarkers()),
_expr = dt * mu0 * sigma * inner(id(A),trans(gradt(V))) );

F0 += integrate(_range=markedelements(mesh, material.meshMarkers()),
_expr = mu0 * sigma * inner(id(A) , idv(Aold)));

M11 += integrate(_range=markedelements(cond_mesh, material.meshMarkers()),

M10 += integrate(_range=markedelements(cond_mesh, material.meshMarkers()),
_expr = mu0 * sigma * inner(idt(A), trans(grad(V))) );

F1 += integrate(_range=markedelements(cond_mesh, material.meshMarkers()),
_expr = mu0 * sigma * inner(idv(Aold), trans(grad(V))) );
}
auto itField = M_modelProps->boundaryConditions().find( "magnetic-potential");
if ( itField != M_modelProps->boundaryConditions().end() )
{
auto mapField = (*itField).second;
auto itType = mapField.find( "Dirichlet" );
if ( itType != mapField.end() )
{
for ( auto const& exAtMarker : (*itType).second )
{
std::string marker = exAtMarker.marker();
auto g = expr<3,1>(exAtMarker.expression());
g.setParameterValues({{"t", t}});
M00 += on(_range=markedfaces(mesh,marker), _rhs=F, _element=*A, _expr= g);
}
}
itType = mapField.find( "DirichletX" );
if ( itType != mapField.end() )
{
for ( auto const& exAtMarker : (*itType).second )
{
std::string marker = exAtMarker.marker();
auto g = expr(exAtMarker.expression());
g.setParameterValues({{"t", t}});
M00 += on(_range=markedfaces(mesh,marker), _rhs=F, _element=(*A)[Component::X], _expr= g);
}
}
itType = mapField.find( "DirichletY" );
if ( itType != mapField.end() )
{
for ( auto const& exAtMarker : (*itType).second )
{
std::string marker = exAtMarker.marker();
auto g = expr(exAtMarker.expression());
g.setParameterValues({{"t", t}});
M00 += on(_range=markedfaces(mesh,marker), _rhs=F, _element=(*A)[Component::Y], _expr= g);
}
}
itType = mapField.find( "DirichletZ" );
if ( itType != mapField.end() )
{
for ( auto const& exAtMarker : (*itType).second )
{
std::string marker = exAtMarker.marker();
auto g = expr(exAtMarker.expression());
g.setParameterValues({{"t", t}});
M00 += on(_range=markedfaces(mesh,marker), _rhs=F, _element=(*A)[Component::Z], _expr= g);
}
}
}
itField = M_modelProps->boundaryConditions().find( "electric-potential");
if ( itField != M_modelProps->boundaryConditions().end() )
{
auto mapField = (*itField).second;
auto itType = mapField.find( "Dirichlet" );
if ( itType != mapField.end() )
{
for ( auto const& exAtMarker : (*itType).second )
{
std::string marker = exAtMarker.marker();
auto g = expr(exAtMarker.expression());
g.setParameterValues({{"t", t}})
M11 += on(_range=markedfaces(cond_mesh,marker), _rhs=F, _element=*V, _expr= g);
}
}
}
...
}``````

Finally, the system is solved by calculating U containing A and V, and then the following time step is carried out.

``    auto result = mybackend->solve( _matrix=M, _rhs=F, _solution=U, _rebuild=true);``

## 2. Results

### 2.1. A solenoidal magnet

The normalized computed electric potential, current and magnetic field at the Origin are plotted bellow:

We use the expected values of the applied electric potential, current and magnetic field for the transient regime (aka t):

Value Unit

$V$

1

V

$I$

135069

A

$B_z(\mathbf{O})$

0.944

T

Here is what we get on feelpp using the block-form program, after normalization by the above values.

To obtain these results, use the command :

``mpirun -np 1 feelpp_mqs_form --config-file cases/quart-turn/quart-turn.cfg --gmsh.hsize 0.05 --pc-type gasm --ksp-monitor=1``

This will produce a csv file in the feel folder associated with this case.

Below is the comparison between getdp and feelpp result :

### 2.2. 2 solenoidal magnets

The normalized computed electric potential, current and magnetic field at the Origin are plotted bellow:

We use the expected values of the applied electric potential, current and magnetic field for the transient regime (aka t):

Value Unit

$V$

1

V

$I$

135069

A

$B_z(\mathbf{O})$

0.850698279

T

Here is what we get on feelpp using the block-form program, after normalization by the above values.

To obtain these results, use the command :

``mpirun -np 1 feelpp_mqs_form --config-file cases/quart-turn/quart-turn2.cfg --gmsh.hsize 0.05 --pc-type gasm --ksp-monitor=1``

This will produce a csv file in the feel folder associated with this case.

Below is the comparison between getdp and feelpp result :