# Gridap.jl field-circuit coupling for magnetostatics

**URL:** <https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069>\
**Category:** Numerics\
**Tags:** package, fem, gridap\
**Created:** [June 20, 2022, 2:06pm UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069 "2022-06-20T14:06:33Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![gwlagerweij](https://avatars.discourse-cdn.com/v4/letter/g/aeb1de/32.png) [@gwlagerweij](https://discourse.julialang.org/u/gwlagerweij)\
**Post date:** [June 20, 2022, 2:06pm UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069/1 "2022-06-20T14:06:33Z")

</div>

Hi all,

I am trying to solve a field-circuit coupled problem using the `Gridap.jl` package.  
I want to simulate the current distribution inside a high-voltage cable, as well as the associated magnetic field.

This means I have two equations I would like to solve simultaneously:

1. A second order differential equation for the magnetic potential in terms of an imposed current density. This equation includes a term for eddy currents (time-harmonic assumption).  
\nabla\cdot(\nu \nabla A\_z) - j\omega\sigma A\_z = -J\_z
2. A circuit equation which imposes a certain current I\_{coil} on the cable.  
\iint\_{coil} J\_z d\Omega = I\_{coil}

I have successfully solved this problem using a home-brew FEM implementation, which manually constructs the linear system with a loop over the elements (see [this Jupyter notebook](https://github.com/gijswl/ee4375_fem/blob/main/fem2d/magnetostatics/cable_skin2.ipynb), and [solution figures](https://github.com/gijswl/ee4375_fem/tree/main/fem2d/magnetostatics/figures) for 3 different frequencies). The results from this code match quite well with those from COMSOL Multiphysics.

I am wondering how to solve this problem using `Gridap.jl`? The first equation is no problem, as this can be solved by just following the first tutorial. However, I am unable to find a suitable starting point for including equation 2 in the documentation, and would greatly appreciate if someone could point me in the right direction.

---

<div class="post-metadata">

**Author:** ![gwlagerweij](https://avatars.discourse-cdn.com/v4/letter/g/aeb1de/32.png) [@gwlagerweij](https://discourse.julialang.org/u/gwlagerweij)\
**Post date:** [July 13, 2022, 11:34am UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069/2 "2022-07-13T11:34:35Z")

</div>

I have been looking into alternatives, and it seems this can be done quite easily using _constraints_ in `Ferrite.jl` ([Constraints · Ferrite.jl](https://ferrite-fem.github.io/Ferrite.jl/dev/manual/constraints/)). I have yet to try this, but it looks promising.

Does anyone know if something similar is possible in `Gridap.jl`?

---

<div class="post-metadata">

**Author:** ![James\_Gorman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/james_gorman/32/16590_2.png) [@James\_Gorman](https://discourse.julialang.org/u/James_Gorman)\
**Post date:** [July 14, 2022, 1:31am UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069/3 "2022-07-14T01:31:46Z")

</div>

Did you look at tutorial 18 ([18 Topology optimization · Gridap tutorials](https://gridap.github.io/Tutorials/dev/pages/t018_TopOptEMFocus/)) by chance? As part of the optimization, they use a radius to filter data (look at the “Topology Optimization” and “Filter and Threshold” sections).

Do you know of a way to rewrite your problem such that you could turn part of it into an optimization problem?

---

<div class="post-metadata">

**Author:** ![gwlagerweij](https://avatars.discourse-cdn.com/v4/letter/g/aeb1de/32.png) [@gwlagerweij](https://discourse.julialang.org/u/gwlagerweij)\
**Post date:** [July 14, 2022, 5:57am UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069/4 "2022-07-14T05:57:00Z")

</div>

Thanks for the pointer! I’ve (very quickly) read through this tutorial and it seems that the section on the [weak form](https://gridap.github.io/Tutorials/dev/pages/t018_TopOptEMFocus/#Weak-form-2) might help in adding an additional equation. I have yet to try this, but will make a post if it succeeds.

> we can then make things simpler by dividing the weak form to a base integral that contains the whole computation cell and an additional integral on the design domain.

As for turning this into an optimization problem, that seems like a waste of computation power: the additional constraint for current stated above _should_ translate to a single extra equation in the linear system, with a marginal effect on the solving of the system.  
Even though the tutorial 18 is labeled “Topology optimization”, I think with the section referred to above, there is no need to actually do any optimization.

---

<div class="post-metadata">

**Author:** ![lijas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lijas/32/20023_2.png) [@lijas](https://discourse.julialang.org/u/lijas)\
**Post date:** [July 14, 2022, 2:40pm UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069/5 "2022-07-14T14:40:45Z")

</div>

> [@gwlagerweij](#):
>
> I have been looking into alternatives, and it seems this can be done quite easily using _constraints_ in `Ferrite.jl` ([Constraints · Ferrite.jl](https://ferrite-fem.github.io/Ferrite.jl/dev/manual/constraints/)). I have yet to try this, but it looks promising.

If you want to use Ferrite, your code in the jupyter notebook should be quite easily translated to Ferrite with the help of the [first tutorial](https://ferrite-fem.github.io/Ferrite.jl/stable/examples/heat_equation/) in the docs. But If I understand you correctly, you also want to implement your second equation using periodic boundary constraints, instead of implementing it in a weak form? It is possible with AffineConstraints, but the tricky part might be to find which dofs to couple on the surfaces…

---

<div class="post-metadata">

**Author:** ![gwlagerweij](https://avatars.discourse-cdn.com/v4/letter/g/aeb1de/32.png) [@gwlagerweij](https://discourse.julialang.org/u/gwlagerweij)\
**Post date:** [July 15, 2022, 8:06am UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069/6 "2022-07-15T08:06:20Z")

</div>

> [@lijas](#):
>
> But If I understand you correctly, you also want to implement your second equation using periodic boundary constraints, instead of implementing it in a weak form? It is possible with AffineConstraints

I agree that the first equation (for the magnetic potential) should be pretty easy to implement using Ferrite. But I’m not sure if I understand what you mean by implementing the second one with periodic boundary constraints.  
I have no preference for the specific implementation, as long as I can do the following:

I have a physical problem: namely, to constrain the total current through a wire to a certain value (I\_{coil}).  
Mathematically, this translates to the second equation stated in my original post.

What do you suggest, implementation-wise, to attain this goal?  
In my own implementation, I have added an additional equation and degree of freedom (additional row and column in the A matrix) for this “constraint”. This is, as far as I am aware, an accepted way of including electrical circuit equations and coupling to the degrees of freedom on the mesh.

---

<div class="post-metadata">

**Author:** ![lijas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lijas/32/20023_2.png) [@lijas](https://discourse.julialang.org/u/lijas)\
**Post date:** [July 18, 2022, 1:02pm UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069/7 "2022-07-18T13:02:41Z")

</div>

I thought you wanted to use periodic boundary conditions because you linked to the it in the Ferrite docs, but perhaps I misunderstood.

Unfortunately I dont know the best method to implement the constraint in Eq2, but the way you have done it with a sort of Lagrange-multiplier methods is pretty standard I think.

---

<div class="post-metadata">

**Author:** ![gwlagerweij](https://avatars.discourse-cdn.com/v4/letter/g/aeb1de/32.png) [@gwlagerweij](https://discourse.julialang.org/u/gwlagerweij)\
**Post date:** [October 15, 2022, 3:07pm UTC](https://discourse.julialang.org/t/gridap-jl-field-circuit-coupling-for-magnetostatics/83069/8 "2022-10-15T15:07:17Z")

</div>

I have managed to sort of implement this using Gridap’s low-level API (as described in [tutorial 13](https://gridap.github.io/Tutorials/dev/pages/t013_poisson_dev_fe)).  
The implementation is available here: [ee4375\_fem\_ta/High-Voltage Cable - Gridap.ipynb at main · gijswl/ee4375\_fem\_ta · GitHub](https://github.com/gijswl/ee4375_fem_ta/blob/main/assignment2_fem2d/High-Voltage%20Cable%20-%20Gridap.ipynb)

After defining the linear and bilinear terms of the weak form, we can use the low-level API to derive the cell contributions and assemble the linear system. This linear system can then be freely modified (e.g., extended with additional rows and columns). The cell contributions to these extra vectors can also be calculated with the low-level API in quite a nice way.