The main subject of this post is to clarify the practical implementation of the Flux Solver in acoustics (might be relevant in other physics as well).
My attempt to compile the latest Elmer source code on Linux Mint 17 64bit did not succeed (based on the error messages, it looks like there are bugs in the source code and/or CMakeLists.txt file. See topic:
viewtopic.php?f=2&t=4922&sid=3410d9d861 ... 324#p17601). Therefore Linux Mint was deleted, and Ubuntu 16.04 64bit was installed in its place. Elmer 8.3 was installed on Ubuntu using the install package from Launchpad (
https://launchpad.net/~elmer-csc-ubuntu ... /+packages). This works, and now I am using the latest Elmer 8.3.
Let’s approach the subject of implementing the FluxSolver from a beginner’s point of view (advanced users don’t need explanations anyway). The initial rough project setup:
- - start a new Elmer project in ElmerGUI
- - import a ready mesh, clean it up, and regenerate it in ElmerGUI
- - set up the model using the menu items in ElmerGUI; but we can not add the FluxSolver this way
- - generate the sif file, save it, and also save the whole project
- - run the solver and check the results in ElmerVTK or Paraview, save the project again.
At this point we have two variable sets in the results, representing the real and imaginary parts of the complex pressure (Pressure wave 1, and Pressure wave 2) in frequency domain.
Our next goal is to compute the acoustic velocity field in ElmerSolver using the FluxSolver. Theoretically it is possible to do this in Paraview, but practically we get wrong results at the in/out boundaries (bug in GradientOfUnstructuredDataSet filter). Therefore, our best option is to compute the velocity using ElmerSolver.
1) After reading the chapter “Flux Computation” in the ElmerModelsManual.pdf, our first attempt is to add the FluxSolver implementation example (given at the end of chapter) to the sif file, and adapt it to computing the acoustic velocity field. Here is the original example:
Code: Select all
Solver 3
Exec Solver = after all
Equation = "flux compute"
Procedure = "FluxSolver" "FluxSolver"
Calculate Flux = Logical True
Flux Variable = String Temperature
Flux Coefficient = String "Heat Conductivity"
Linear System Solver = "Iterative"
Linear System Iterative Method = "cg"
Linear System Preconditioning = ILU0
Linear System Residual Output = 10
Linear System Max Iterations = Integer 500
Linear System Convergence Tolerance = 1.0e-10
End
The line “Flux Variable = String Temperature” is confusing, because the flux variable supposed to be a vector, while the Temperature is a scalar. Performing a search in the pdf for “Flux Variable” we get 5 hits. The one on page 198 defines “Flux Variable” as:
“Flux Variable String
This gives the name of the flux variable used to compute the source term. Note that this must be the name of a vector field such as Velocity.”
Our logical thought is that this instruction supposed to define the name of the FluxSolver’s output variable, which in our case is the velocity vector. Let’s name it “AcousticVelocity”. So the discussed line will be modified to look like this:
Flux Variable = String “AcousticVelocity”
2) At this point we assume that the FluxSolver will take the main solver output (which is the pressure field), calculate it’s gradient, and multiply this with the Flux Coefficient of "Heat Conductivity".
In our case the coefficient supposed to be 1/(Density*Angular Frequency) which is 3.85e-4 here.
Let’s make it simple now, and just assign a value of 3.85e-4 to the “Heat Conductivity” in the section “Material 1” of sif file. We can do this as long as the real value of “Heat Conductivity” is not used in our project.
Heat Conductivity = 3.85e-4
Next, let’s enable the FluxSolver in the “Equation 1” section:
Code: Select all
Equation 1
Name = "Helmholtz"
Angular Frequency = 2155
Active Solvers(3) = 2 3 1
End
At present Solver 2 is the main “HelmholtzSolver”, Solver 3 is the “FluxSolver”, and Solver 1 is the "ResultOutputSolver". Just to prevent unnecessary problems, the sequential order of solvers 2 3 1 is given to execute the main solver first, then the FluxSolver, and finally the ResultOutputSolver.
Why not renumber the solvers to make the sequence 1 2 3 do the same? Because the ElmerGUI originally has generated an inverted sequential order, making the main HelmholtzSolver to be #2 and the ResultOutputSolver to be #1, while using “Active Solvers(2) = 2 1” sequence. Therefore, let’s keep this unchanged for now, and just insert our FluxSolver as #3 between the two original solvers, thus getting the sequence 2 3 1.
But, for clarity: according to the manual the sequential order of solvers is unimportant as long as the values of the “Exec Solver =” variable are correctly set in every solver section.
At this point the sif file looks like this:
Code: Select all
Header
CHECK KEYWORDS Warn
Mesh DB "." "."
Include Path ""
Results Directory ""
End
Simulation
Max Output Level = 10
Coordinate System = Cartesian
Coordinate Mapping(3) = 1 2 3
Simulation Type = Steady state
Steady State Max Iterations = 1
Output Intervals = 1
Timestepping Method = BDF
BDF Order = 1
Solver Input File = tube.sif
Post File = tube.ep
End
Constants
Gravity(4) = 0 -1 0 9.82
Stefan Boltzmann = 5.67e-08
Permittivity of Vacuum = 8.8542e-12
Boltzmann Constant = 1.3807e-23
Unit Charge = 1.602e-19
End
Body 1
Target Bodies(1) = 1
Name = "Body Property 1"
Equation = 1
Material = 1
Initial condition = 1
End
Solver 2
Equation = Helmholtz Equation
Procedure = "HelmholtzSolve" "HelmholtzSolver"
Variable = -dofs 2 Pressure Wave
Exec Solver = Always
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-8
Nonlinear System Max Iterations = 20
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Direct
Linear System Direct Method = Banded
End
Solver 3
Exec Solver = after all
Equation = "flux compute"
Procedure = "FluxSolver" "FluxSolver"
Calculate Flux = Logical True
!Flux Variable = String Temperature
Flux Variable = String “AcousticVelocity”
Flux Coefficient = String "Heat Conductivity"
Linear System Solver = "Iterative"
Linear System Iterative Method = "cg"
Linear System Preconditioning = ILU0
Linear System Residual Output = 10
Linear System Max Iterations = Integer 500
Linear System Convergence Tolerance = 1.0e-10
End
Solver 1
Equation = Result Output
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Save Geometry Ids = False
Output File Name = tube
Output Format = Vtu
Scalar Field 2 = Pressure wave 2
Scalar Field 1 = Pressure wave 1
Exec Solver = After Simulation
End
Equation 1
Name = "Helmholtz"
Angular Frequency = 2155
Active Solvers(3) = 2 3 1
End
Material 1
Name = "Air (room temperature)"
Viscosity = 1.983e-5
Heat expansion Coefficient = 3.43e-3
Heat Conductivity = 3.85e-4
Relative Permittivity = 1.00059
Sound speed = 343.0
Heat Capacity = 1005.0
Density = 1.205
End
Initial Condition 1
Name = "InitialCondition 1"
Pressure Wave 1 = 0
Pressure Wave 2 = 0
End
Boundary Condition 1
Target Boundaries(1) = 6
Name = "In"
Pressure Wave 2 = 0
Pressure Wave 1 = 100
End
Boundary Condition 2
Target Boundaries(4) = 2 3 4 5
Name = "Walls"
Wave Flux 1 = 0
Wave Flux 2 = 0
End
Boundary Condition 3
Target Boundaries(1) = 1
Name = "Out"
Wave impedance 1 = 343
Wave impedance 2 = 0
End
Now we save the modified sif file as backup with a name “tube_step1.sif0”, and also save it as “tube.sif”, which will be used by the solver. After running the solver, and opening the solution in ElmerVTK postprocessor (flux was not exported yet into the results file for Paraview), 3 additional output variable sets appear in the solution, all of them empty (trivially zero):
acousticvelocity.flux_x
acousticvelocity.flux_y
acousticvelocity.flux_z
This is in agreement with the relevant report lines in the Solver log:
“SolveLinearSystem: Assuming real valued linear system
SolveSystem: Solution trivially zero!”
3) Apparently the FluxSolver did not use the main HelmhotlzSolver’s output variables of “Pressure Wave” for the flux computation as we have assumed it should have done (based on the example code, which did not specify the input variable explicitly). OK, let's specify this now, and add this line to the section “solver 3”:
Target Variable = String "Pressure Wave 2"
and modify the name of the output velocity (flux) name (the real part of the complex velocity is calculated from the gradient of the imaginary part of the complex pressure; see the equations given earlier):
Flux Variable = String “AcousticVelocity1”
Save the modifications in the sif file, and rerun the solver. The flux is still trivially zero.
4) What is wrong? A hunch tells me “let's remove the “Flux Variable” instruction, because it was not on the list of possible keywords for FluxSolver in the manual”. Therefore remove the:
Flux Variable = String “AcousticVelocity1”
from the sif file, save it, and run the solver again.
5) Now the log file shows that the ResultOutputSolver is executed before the FluxSolver, but that should not be the cause of an empty flux variable in the .ep file. However, the “after all” type of FluxSolver execution makes no sense, because it has to be executed before saving. Let's try to fix the execution sequence with these modifications (chosen from the list of possible options given in ElmerSolverManual.pdf: never, always, before all, before timestep, after timestep, before saving, after saving, after all), which would make sense:
In HelmholtzSolver section let's keep the original:
Exec Solver = String "always"
In FluxSolver section:
Exec Solver = String "before saving"
In ResultOutputSolver section:
Exec Solver = String "before saving"
(By the way, the “After Simulation” option generated by the ElmerGUI is not present in the ElmerSolverManual.pdf)
Save and run. Now the solver log shows that the FluxSolver started the computation, but it stopped with an error message:
“NUMERICAL ERROR:: IterSolve: Failed convergence tolerances.
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL”
When starting ElmerVTK post processor to view the results, ElmerGUI freezes. Force quit the application, restart ElmerGUI, load the project, open the last major sif file backup in text editor, and save it as tube.sif. Now we can continue the debugging process.
6) Decreased the convergence tolerance from 1e-10 to:
Linear System Convergence Tolerance = 1.0e-5
but this did not fix the problem. Then wasted some more hair by trying different settings of the iterative solver, but I don’t want to bore you with all that here.
7) At this stage our aim is not the experimentation with solver settings, but rather just getting the flux solver to work with any solver. Let’s use therefore the simplest direct solver that worked nicely in the HelmholtzSolver. Replace the “Linear System...” lines in FluxSolver section with these:
Linear System Solver = Direct
Linear System Direct Method = Banded
Save, and run the solver. Success! Finally the FluxSolver gave correct results for the real components of acoustic velocity field. Now we have got these solver output variables in the tube.ep file with correct values:
Pressure wave 1
Pressure wave 2
Pressure wave 2 flux_x
Pressure wave 2 flux_y
Pressure wave 2 flux_z
8) Next step is to get the flux exported into the tube.vtu file for Paraview. Let's add this line to the ResultOutputSolver section:
Vector Field 1 = Pressure wave 2 flux
Save and run. If we open the tube0001.vtu output file in Paraview, the variable “Pressure wave 2 flux” is present, but it is trivially zero.
9) We see in the solver log that this happens because the vtu file is written before the flux is calculated. To avoid further unnecessary confusion, let's remove all the “Exec Solver = ...” lines from the sif file, and renumber the sequential order of the solvers in the order we want them to be executed:
Solver 1 → HelmholtzSolver
Solver 2 → FluxSolver
Solver 3 → ResultOutputSolver
Active Solvers(3) = 1 2 3
Save the sif file, run the solver. Now the “Pressure wave 2 flux” is not only present in Paraview, but it also contains correct values for all 3 components. In hindsight, we should have avoided using the confusing “Exec Solver = ...” instructions from the very beginning, because a strictly defined sequence by numbering is a clearer and simpler way of controlling what should be computed first, and what later.
10) Next, create one more FluxSolver section for the imaginary part of the velocity:
Code: Select all
Solver 3
Equation = "flux compute 2"
Procedure = "FluxSolver" "FluxSolver"
Calculate Flux = Logical True
Target Variable = String "Pressure Wave 1"
Flux Coefficient = String "Heat Conductivity"
Linear System Solver = Direct
Linear System Direct Method = Banded
End
Also implement these modifications:
- In the first FluxSolver change its name to:
- Equation = "flux compute 1"
- renumber the ResultOutputSolver section to Solver 4, and add this line to it:
- Vector Field 2 = Pressure wave 1 flux
- Active Solvers(4) = 1 2 3 4
Save the sif file, run the solver. The solution now contains both Re and Im components of pressure scalar field, and also of the velocity vector field, with correct values in both output files for VTKpost and Paraview as well. The sif file at this point looks like this:
Code: Select all
Header
CHECK KEYWORDS Warn
Mesh DB "." "."
Include Path ""
Results Directory ""
End
Simulation
Max Output Level = 10
Coordinate System = Cartesian
Coordinate Mapping(3) = 1 2 3
Simulation Type = Steady state
Steady State Max Iterations = 1
Output Intervals = 1
Timestepping Method = BDF
BDF Order = 1
Solver Input File = tube.sif
Post File = tube.ep
End
Constants
Gravity(4) = 0 -1 0 9.82
Stefan Boltzmann = 5.67e-08
Permittivity of Vacuum = 8.8542e-12
Boltzmann Constant = 1.3807e-23
Unit Charge = 1.602e-19
End
Body 1
Target Bodies(1) = 1
Name = "Body Property 1"
Equation = 1
Material = 1
Initial condition = 1
End
Solver 1
Equation = Helmholtz Equation
Procedure = "HelmholtzSolve" "HelmholtzSolver"
Variable = -dofs 2 Pressure Wave
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-8
Nonlinear System Max Iterations = 20
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Direct
Linear System Direct Method = Banded
End
Solver 2
Equation = "flux compute 1"
Procedure = "FluxSolver" "FluxSolver"
Calculate Flux = Logical True
Target Variable = String "Pressure Wave 2"
Flux Coefficient = String "Heat Conductivity"
Linear System Solver = Direct
Linear System Direct Method = Banded
End
Solver 3
Equation = "flux compute 2"
Procedure = "FluxSolver" "FluxSolver"
Calculate Flux = Logical True
Target Variable = String "Pressure Wave 1"
Flux Coefficient = String "Heat Conductivity"
Linear System Solver = Direct
Linear System Direct Method = Banded
End
Solver 4
Equation = Result Output
Procedure = "ResultOutputSolve" "ResultOutputSolver"
Save Geometry Ids = False
Output File Name = tube
Output Format = Vtu
Scalar Field 1 = Pressure wave 1
Scalar Field 2 = Pressure wave 2
Vector Field 1 = Pressure wave 2 flux
Vector Field 2 = Pressure wave 1 flux
End
Equation 1
Name = "Helmholtz"
Angular Frequency = 2155
Active Solvers(4) = 1 2 3 4
End
Material 1
Name = "Air (room temperature)"
Viscosity = 1.983e-5
Heat expansion Coefficient = 3.43e-3
Heat Conductivity = 3.85e-4
Relative Permittivity = 1.00059
Sound speed = 343.0
Heat Capacity = 1005.0
Density = 1.205
End
Initial Condition 1
Name = "InitialCondition 1"
Pressure Wave 1 = 0
Pressure Wave 2 = 0
End
Boundary Condition 1
Target Boundaries(1) = 6
Name = "In"
Pressure Wave 2 = 0
Pressure Wave 1 = 100
End
Boundary Condition 2
Target Boundaries(4) = 2 3 4 5
Name = "Walls"
Wave Flux 1 = 0
Wave Flux 2 = 0
End
Boundary Condition 3
Target Boundaries(1) = 1
Name = "Out"
Wave impedance 1 = 343
Wave impedance 2 = 0
End
11) Next step is to find at least one iterative solver setting that could be used in the FluxSolver instead of the direct solver. Here is a version, which is based on a code snippet of Matthias from another thread.
Delete these two lines from the FluxSolver sections:
Linear System Solver = Direct
Linear System Direct Method = Banded
and replace them with these:
Code: Select all
Stabilize = True
Bubbles = False
Lumped Mass Matrix = False
Optimize Bandwidth = True
Steady State Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-7
Nonlinear System Max Iterations = 20
Nonlinear System Newton After Iterations = 3
Nonlinear System Newton After Tolerance = 1.0e-3
Nonlinear System Relaxation Factor = 1
Linear System Solver = Iterative
Linear System Iterative Method = BiCGStab
Linear System Max Iterations = 500
Linear System Convergence Tolerance = 1.0e-10
BiCGstabl polynomial degree = 2
Linear System Preconditioning = Diagonal
Linear System ILUT Tolerance = 1.0e-3
Linear System Abort Not Converged = False
Linear System Residual Output = 1
Linear System Precondition Recompute = 1
Save the sif and run the solver. Now the results still look good, but in the solver log there are two of these warnings:
“88 0.3180E-38
NUMERICAL WARNING:: IterSolve: Failed convergence tolerances.”
12) Experimented with lower convergence tolerance values (implementing only one change at a time) like:
Linear System Convergence Tolerance = 1.0e-5
Nonlinear System Convergence Tolerance = 1.0e-5
Steady State Convergence Tolerance = 1.0e-3
with no luck, the warnings were still there. Therefore, at this point I am giving up the fiddling with the iterative solver, in hope that someone who really understands what all these parameters actually do, could show some settings that work perfectly.
Summary:
- - The “Flux Variable = ...” must not be present in the FluxSolver section, because then the solver breaks down, giving no flux values (this is a mistake in the manual). With other words we can't define what the name of FluxSolver’s output variable should be.
- - At the same time, we must use “Target Variable = ...” to define the input scalar field for the FluxSolver (which was not present in the manual’s example).
- - The simplest and least confusing is not to use “Exec Solver = ...” instructions, but use sequentially correct solver numbering and order instead, to make sure that the flux is calculated after the main solver, and before the vtu file gets written.
- - If there are more than one FluxSolver sections, then each one must have a unique name in “Equation = ...” instructions, otherwise the solver won't work.
- - For the time being the setting “Flux Coefficient = String "Heat Conductivity"” works, but this is not an ideal solution. Does anybody know how to use a different freely chosen parameter (or variable) name, instead of Heat Conductivity? The best would be if we would not have to calculate this value of 3.85e-4 manually, but let the solver calculate it from this equation: 1/(Density*Angular Frequency).
- - The iterative solver given in the example of the manual does not work. The direct solver does. The BiCGStab works as well, but with incomplete convergence. Looking for suggestions about which iterative solvers would be ideal for this purpose, with what settings.
Thanks for your patience if you have read this long post. The conclusions and code samples aim to help beginners. The description of the lengthy process of “figuring it out” will hopefully make the point to the experts and developers, how much collective time and nerve can be unnecessarily wasted on rediscovering simple things, due to the lack of detailed, clear, and accurate instructions for beginner users.
If somebody discovers (or a priory knows) how to use a feature, which has not been accurately described in detail in manuals (or elsewhere), please sacrifice some of your time to write it down, and share it with everybody! There is a huge difference between one person wasting a day or two on figuring out a feature, and hundreds of users doing the same, wasting thousands of collective hours on rediscovering the wheel. Such explanations supposed to be also collected into a single document for easy access, instead of leaving them scattered here and there, making it less likely for people in need to find them, when they need them.
Looking forward to your comments and solutions.
Regards,
Zoltan