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
Edge degrees of freedom
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 curlconforming fields simultaneously hasn't been exemplified yet, although I have been having a longlasting 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 twodimensional 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
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 curlconforming fields simultaneously hasn't been exemplified yet, although I have been having a longlasting 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 twodimensional 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
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
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
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 complexvalued 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
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
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
Variable DOFs = 4
Element = "n:0 e:1"
should in principle serve as well, with a chance of getting more useful aliases.
 Mika
Re: Edge degrees of freedom
Thanks for the help Mika!
I was able to implement the weak ATformulation 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 nonconductive (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 TDOFs in the air.
As far as I can understand, in Elmer it is not possible to define a bodywise 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 ATSolver module or would such an option require changes within the Elmer source code?
I made a comparison of ATformulation and WhitneyAVharmonicSolver which I previously adjusted so that the AVformulation 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 ATformulation compared to AVformulation.
* 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 TDOFs, 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 nonzero values, although it should contain only nonzero 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 nonzero 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
I was able to implement the weak ATformulation 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 nonconductive (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 TDOFs in the air.
As far as I can understand, in Elmer it is not possible to define a bodywise 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 ATSolver module or would such an option require changes within the Elmer source code?
I made a comparison of ATformulation and WhitneyAVharmonicSolver which I previously adjusted so that the AVformulation 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 ATformulation compared to AVformulation.
* 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 TDOFs, 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 nonzero values, although it should contain only nonzero 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 nonzero 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
Re: Edge degrees of freedom
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.
The cost of a matrixvector 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 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?
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 sparsematrix representation.
Best regards,
Mika

 Site Admin
 Posts: 4133
 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:
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 matrixvector 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
Great that you're trying out new formulations!
There is a keyword:
Code: Select all
Linear System Remove Zeros = Logical True
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 matrixvector 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
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 AVformulation and the ATformulation by specifying a completely nonconductive domain (air only). In that case, the two formulations should be equivalent. What I get is that the STIFF matrix in the ATformulation has twice as many elements as the STIFF matrix in the AVformulation, but both formulations have STIFF matrices with the same number of nonzero 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 ATformulation, one iterative step takes almost twice as long as in the case of the AVformulation. It seems to me that the reason for twice the duration of the ATsimulation 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 nonzero elements, more precisely about 2,600,000 nonzero elements. I am fairly convinced that removing rows related to Tdegrees 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 AVformulation variable has 4 degrees of freedom, while the AVformulation 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
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 AVformulation and the ATformulation by specifying a completely nonconductive domain (air only). In that case, the two formulations should be equivalent. What I get is that the STIFF matrix in the ATformulation has twice as many elements as the STIFF matrix in the AVformulation, but both formulations have STIFF matrices with the same number of nonzero 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 ATformulation, one iterative step takes almost twice as long as in the case of the AVformulation. It seems to me that the reason for twice the duration of the ATsimulation 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 nonzero elements, more precisely about 2,600,000 nonzero elements. I am fairly convinced that removing rows related to Tdegrees 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 AVformulation variable has 4 degrees of freedom, while the AVformulation 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