## Edge degrees of freedom

Discussion about coding and new developments
spanda
Posts: 11
Joined: 19 Feb 2018, 14:24
Antispam: Yes

### Edge degrees of freedom

Hi

Is it possible to use a variable that has more than one degrees of freedom per edge? By this I do not mean edge basic functions of the second degree.

For example, is it possible to use a weak formulation containing E and H vector fields (where E represents the electric field vector and H represents the magnetic field vector) such that both vector fields are approximated by edge elements.

I guess it would then look like this in a SIF file:
----------------------
Solver 1
..
..
Variable = Field[ E re:1 E im:1 H re:1 H im:1] !...harmonic case
Variable DOFs= 4
Element = "n:0 e:1"
..
..
End
-------------------------------

I am trying to develop similar formulation in the 3D case, but the solver does not seem to know that the variable has 4 edge degrees of freedom, i.e. real and imag part for both E and H.

For example:
When I print final solution for an element that is a tetrahedron (so it has 6 edges), by using the GetLocalSolution() subroutine, the solutuion vector SOL has only 14 values ​​instead of 24 values ​​(4 EdgeDofs for each of the 6 edges of the tetrahedron).

It seems to me that it offers me real an imaginary part for each of the 6 edges plus 2 strange values (maybe some unwanted nodal Dofs), from which it follows:
2 * 6Edges + 2dofs = 14 values ​​in SOL.

Thanks and best regards
--Spanda
mika
Posts: 202
Joined: 15 Sep 2009, 07:44

### Re: Edge degrees of freedom

Hi,

Elmer should offer pretty much general support for associating degrees of freedom with different geometric entities of a background mesh, so this should be possible to do. Nevertheless the case of approximating two curl-conforming fields simultaneously hasn't been exemplified yet, although I have been having a long-lasting plan to do something similar.

How many DOF values the subroutine GetLocalSolution(X) returns should depend on the return value of the function GetElementDOFs, if the size of the result array X(:,:) doesn't make a restriction. I'd first check that the return value of the function GetElementDOFs corresponds to your specification of DOFs and also check that the two-dimensional array X has sufficient dimensions.

I made some simple experiments with element definitions which could be utilized in this case, but I didn't see anything suspicious outright.

-- Mika
spanda
Posts: 11
Joined: 19 Feb 2018, 14:24
Antispam: Yes

### Re: Edge degrees of freedom

Thanks for the reply Mika!

To be more precise, this is the AT weak formulation shown in the paper by Albertz and Henneberger "Calculation of 3D Eddy Current Fields Using Both Electric and Magnetic Vector Potential in Conducting Regions", IEEE Trans. Magn.
In this paper, both A and T potentials are approximated using edge DOFs.

Do you have any advice on how to properly give instructions to Elmer so that it works with two edgeDOFs on each edge, both with real and imaginary part?

After examining the GetElementDOFs() function, it seems to me that EDOFs = Solver% Mesh% Edges (Element% EdgeIndexes (j))% BDOFs parameter should be equal to 2, or alternatively parameter EdgeDOFs = Solver % Mesh % MaxEdgeDOFs should be equal to 2.
I suppose then the following should be specified in the SIF file:
------------
Solver 1
..
..
Variable = Field[ E re:1 E im:1 H re:1 H im:1] !...harmonic case
Variable DOFs= 2
Element = "n:0 e:2"
..
..
End
--------------------
With this approach parameters EDOFs = Solver% Mesh% Edges (Element% EdgeIndexes (j))% BDOFs and EdgeDOFs = Solver % Mesh % MaxEdgeDOFs are set to 2.

So, Element = "n:0 e:2" command will specify two edgeDOFs for each edge, and Variable DOFs= 2 command is there due to complex values for the DefaultUpdateEquationsC subroutine to correctly assemble the STIFF matrix and the RHS vector.

But, my simulation doesn't converge so maaybe this approach is incorrect

Am I wrong? Is it enough to use only a SIF file as described above or it is necessary to add some unmentioned (and unknown to me) commands in the NewSolver_init (initialization phase) in order to properly set the number of edgeDOFs to 2?

Thank you Mika for the help

--Spanda
mika
Posts: 202
Joined: 15 Sep 2009, 07:44

### Re: Edge degrees of freedom

I'd also set "Variable DOFs = 2" so that the two components would represent the real and imaginary parts. This is the way how Elmer has traditionally handled the approximation of complex-valued fields and following this convention offers a better way to utilize general library routines. The element definition Element = "n:0 e:2" should also give the right set of degrees of freedom for tetrahedra.

Using aliases in the variable name definition needs however some care. Now it will be difficult to get an alias which would follow the splitting of the unknown to A and T. When one uses "Variable DOFs = 2", aliases could be defined by setting for example

Variable = AT[AT re:1 AT im:1]

Then "AT re" would refer to a variable which contains values for representing the real part of both A and T. Unfortunately this is rather abstract. A postprocessing routine could be written to create "A Re", "A Im", "T Re" and "T Im" from "AT". A tailored subroutine for setting BCs seems also necessary due to abstract nature of the variable.

When two DOFs per real/imaginary component are defined on an edge, the first DOF could be used for representing A and the second for representing T. The convention may also be something else, but the convention must be followed consistently when writing code.

Best regards,
Mika
mika
Posts: 202
Joined: 15 Sep 2009, 07:44

### Re: Edge degrees of freedom

I add that the way to get the right number of DOFs is not unique nevertheless. I think the definitions

Variable DOFs = 4
Element = "n:0 e:1"

should in principle serve as well, with a chance of getting more useful aliases.

-- Mika
spanda
Posts: 11
Joined: 19 Feb 2018, 14:24
Antispam: Yes

### Re: Edge degrees of freedom

Thanks for the help Mika!

I was able to implement the weak AT-formulation shown in the paper by Albertz and Henneberger "Calculation of 3D Eddy Current Fields Using Both Electric and Magnetic Vector Potential in Conducting Regions", IEEE Trans. Magn.

I used the following approach:

Variable = P [A re: 1 A im: 1 T re: 1 T im: 1]
Variable DOFs = 4
Element = "n: 0 e: 1"

Elmer successfully computes both A and T vector potentials approximated with edge DOFs. I should have reported it here a year ago.

Now, I have a few questions. I apologize for the length of the text.
In the case where my problem domain consists of both conductive (iron) and non-conductive (air) regions, there is of course a problem with redundant degrees of freedom in the air. Since the air is an eddy current free region, I only need the magnetic vector potential A, i.e. the DOFs for T are redundant and the simulation takes almost twice as long. Therefore, I would like to eliminate T-DOFs in the air.
As far as I can understand, in Elmer it is not possible to define a body-wise variable, i.e., to specify something like:
.
.
Body 1
Equation = 1
Material = 1! Iron
End

Body 2
Equation = 2
Material = 2! Air
End

Equation 1
Active Solvers (2) = 1 2
Variable = P [A re: 1 A im: 1 T re: 1 T im: 1]
End

Equation 2
Active Solvers (2) = 1 2
Variable = A [A re: 1 A im: 1]
End

.
.
because in general ElmerSolver can use only one variable within the whole domain? Would it be possible to implement something similar within the AT-Solver module or would such an option require changes within the Elmer source code?

I made a comparison of AT-formulation and WhitneyAVharmonicSolver which I previously adjusted so that the AV-formulation is symmetrical so that I can use the CG algorithm in both cases. I got the following results:
* One CG iteration step takes almost twice as long in the case of AT-formulation compared to AV-formulation.
* Number of total degrees of freedom for AT = 308 668,...I got it by printing SIZE (ATSolver% Variable% Values).
* Number of total degrees of freedom for AV = 181 642,...I got it by printing SIZE (AVSolver% Variable% Values).
* In both cases the CRS matrix format is used. The size of the Solver%Matrix%Values ​​vector is twice as large, ie:
..... for AT I got ... SIZE(ATSolver%Matrix%Values) = 22 362 096
..... for AV I got ... SIZE(AVSolver%Matrix%Values) = 11 648 524
I assume the size of the ATSolver%Matrix%Values vector and not the ATSolver%Variable%Values vector should actually be the main reason for the longer duration of the one CG simulation step of the ATSolver compared to ATSolver (i.e WhitneyAVHarmonicSolver). Am I right?

So what interests me is, is there perhaps a different strategy for eliminating redundant T-DOFs, in order to shorten the simulation duration. Maybe through manipulating the ATSolver%Matrix%Values ​​vector if it is not already possible to have bodywise defined variables? Do you have any advice?

Some additional notes:
I placed all of the above printing commands in two different positions in the code. The first position is immediately before the CALL DefaultDirichletBCs() command and the second position is immediately before the execution of the Norm = DefaultSolve() command.
In both cases, I got the same numbers, which was a bit unexpected for me, so I printed values ​​from the Solver%Matrix%Values vector to see if DirichletBCs change coefficients in STIFF matrix at all. And, there was a change in the coefficients after CALL DefaultDirichletBCs() command but it did not affect the size of the vector Solver%Matrix%Values (in both AT and WhitneyAV case) becouse Solver%Matrix%Values vector contained zeros in addition to non-zero values, although it should contain only non-zero values ​​as far as I know. Does ElmerSolver only subsequently (within the DefaultSolve() routine) eliminate zeros from the Solver% Matrix%Values vector or is there something wrong with the CRS matrix procedure?
I then tried to count non-zero values ​​in the Solver% Matrix% Values vector using the Count(Solver%Matrix%Values) /= 0.0 command and got 4 637 752 instead of 11 648 524. Is this somehow expected? By the way, simulation results seems fine. Also, size of Matrix%Cols vector was equal to size of Matrix%Values vector.

Thanks and best regards
--Spanda
mika
Posts: 202
Joined: 15 Sep 2009, 07:44

### Re: Edge degrees of freedom

spanda wrote: 03 Jun 2021, 21:49 in general ElmerSolver can use only one variable within the whole domain? Would it be possible to implement something similar within the AT-Solver module or would such an option require changes within the Elmer source code?
This is not possible, since the unknown of a PDE model is supposed to have the same components everywhere. Enabling this would require changes to the general utilities and this might be a complicated task. When unnecessary DOFs exist, the system matrix is usually modified such that the corresponding rows have diagonal entries only to obtain the zero solution (the corresponding RHS entry is set to be zero). This strategy is the same as what is applied when adding zero Dirichlet BCs.

I'm not sure how much these superfluous constraint rows affect the solution time, but I believe their effect may not be the most important factor. I expect that a key to efficient solution is the choice/design of a preconditioner.
spanda wrote: 03 Jun 2021, 21:49 I assume the size of the ATSolver%Matrix%Values vector and not the ATSolver%Variable%Values vector should actually be the main reason for the longer duration of the one CG simulation step of the ATSolver compared to ATSolver (i.e WhitneyAVHarmonicSolver). Am I right?
The cost of a matrix-vector product should be proportional to the number of the matrix entries, but the cost of the application of the preconditioner may also be an important factor here.
spanda wrote: 03 Jun 2021, 21:49 I then tried to count non-zero values ​​in the Solver% Matrix% Values vector using the Count(Solver%Matrix%Values) /= 0.0 command and got 4 637 752 instead of 11 648 524. Is this somehow expected?
A local matrix (elementwise contribution) usually contains some zero entries and they become assembled into the global system. The assembled global matrix is thus not a perfect sparse-matrix representation.

Best regards,
Mika
raback
Site Admin
Posts: 4173
Joined: 22 Aug 2009, 11:57
Antispam: Yes
Location: Espoo, Finland
Contact:

### Re: Edge degrees of freedom

Hi

Great that you're trying out new formulations!

There is a keyword:

Code: Select all

``````Linear System Remove Zeros = Logical True
``````
That strips away all zeros. The caveat is that the CRS matrix structure is different when coming back to assembly the matrix next time. But this might give a clue whether there is much to be gained. This might not work with all linear solvers and preconditioners either. Therefore it is not adverticed either.

As Mika stated, also I think there is little to be gained in eliminating some zeros coming from Dirichlet BCs. On the other hand, would there be zero blocks in the matrix equation the case is quite different. The number of entries determine the cost of matrix-vector product. It does not help to test out the zeros, because it costs easily at least as much as just going on with the product.

Are there zero blocks in your system?

-Peter
spanda
Posts: 11
Joined: 19 Feb 2018, 14:24
Antispam: Yes

### Re: Edge degrees of freedom

Hi Peter
Thanks for the help!

On your advice, I tried keyword: Linear System Remove Zeros = Logical True

In WhitneyAVHarmonicSolver I didn't notice any difference whether the keyword is True or False. On the other hand, in WhitneyAVSolver the simulation diverges when the keyword is True.

I see that the keyword ultimately calls the CRS_RemoveZeros() subroutine which is called in the WhitneyAVSolver via the DefaultFinishBulkAssembly() subroutine.
As far as I can see, the DefaultFinishBulkAssembly() subroutine is not present in WhitneyAVHarmonicSolver, nor is the CRS_RemoveZeros() subroutine called in any alternative way. I guess that's why that keyword has no effect in WhitneyAVHarmonicSolver, regardless of whether it's True or False. Is there a reason why the CRS_RemoveZeros() subroutine is omitted from the WhitneyAVHarmonicSolver?
I added the CRS_RemoveZeros() subroutine to my own module but the simulation also diverged. The CRS_RemoveZeros() routine seems to remove some zeros that it probably shouldn't remove. Maybe I should write a subroutine based on the CRS_RemoveZeros() subroutine that would only remove rows that are related to T degrees of freedom in the air.

In my case, there are almost half of the rows in the STIFF matrix that are redundant (rows with only zero elements). I don’t know if such a matrix can be said to have zero blocks?

I made an interesting test. I made a comparison between the AV-formulation and the AT-formulation by specifying a completely non-conductive domain (air only). In that case, the two formulations should be equivalent. What I get is that the STIFF matrix in the AT-formulation has twice as many elements as the STIFF matrix in the AV-formulation, but both formulations have STIFF matrices with the same number of non-zero elements that are also equal, as expected. Also, both simulations converge in the same number of steps (266 iterative steps) and have the same residual norms. Where these two simulations differ is the duration of the individual iterative step. In the case of the AT-formulation, one iterative step takes almost twice as long as in the case of the AV-formulation. It seems to me that the reason for twice the duration of the AT-simulation is precisely the size of the Matrix%Values ​​vector, where the ATSolver%Matrix%Values ​​vector has about 22,000,000 elements, whereas the AVSolver%Matrix%Values ​​vector has about 11,500,000 elements. As I said before, both vectors have the same number of non-zero elements, more precisely about 2,600,000 non-zero elements. I am fairly convinced that removing rows related to T-degrees of freedom (in air) in the STIFF matrix (i.e., in the ATSolver%Matrix%Values ​​vector) would generally result in equal simulation durations.

Note: Although the AV-formulation variable has 4 degrees of freedom, while the AV-formulation variable has 6 degrees of freedom (so 1.5 times more), there are many more edges than nodes in a tetrahedron mesh, so I guess that's why the ATSolver%Matrix%Values ​​vector is twice as large as the AVSolver%Matrix%Values ​​vector.

Best regards
Spanda