Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Numerical methods and mathematical models of Elmer
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by Z.Losonc »

The “Simple Pipe Acoustics” project posted by JoeD motivated me to try using Elmer for acoustic simulations. I have never used Elemer before, so please be patient with me. It would be really great if we could create a collection of sorely needed acoustics tutorials that drill down to the nitty-gritty in a way that beginners could also understand.

I will use a modified version of Joe’s project here, which simulates “the sound propagation within a 1 m long pipe of 0.1 m diameter filled with air, having acoustically hard walls. The input and output openings are matched boundaries to avoid reflections. At the input we have an incident plane wave of p0=100 Pa amplitude.” The process of setting up the simulation is not a trivial task for beginners, and it will need to be described at some point, but let’s focus first on solving a few problems.
tube_p1_mesh.png
tube_p1_mesh.png (52.31 KiB) Viewed 8696 times
The output of the solver is a complex pressure field of which the real part is named “Pressure Wave 1” and the imaginary part is “Pressure Wave 2”. With this we can create nice color images in Paraview and also see the pressure values at any point, but for most practical purposes this is not enough. In most cases we want to know the absolute pressure field, particle velocity field, the radiation intensity field, the power radiated through a surface, the reflected wave amplitude/intensity, SPL field, the far field radiation pattern etc. and also to create some animations. Let’s perform some post processing here, and calculate at least some of the mentioned variables.

1) First of all we would like to create an animation that shows the movement of the pressure field in the tube, so that we could gain a visual insight into the existence of standing waves, propagating waves, the direction of motion etc. In the ElmerModelsManual.pdf, “Helmholtz Solver” chapter the following equation (9.3) shows how to calculate the pressure field at a specific time t using the complex pressure variables provided by the solver:
Eq_1.png
Eq_1.png (3.15 KiB) Viewed 8684 times
The om*t is the phase angle of the complex pressure, so for the purpose of the animation we can simplify this by replacing it with an expression for a specific phase angle fi. This can be easily implemented in Paraview by entering the following formula into the calculator box (here fi=90deg):

p(fi)=pressure wave 1*cos(90*3.14/180)-pressure wave 2*sin(90*3.14/180)

This will create a pressure field that is shifted by 90 degrees in phase. To generate the field at 10 degrees, just replace the number 90 in the above equation and press the button [Apply]. If you save the image after each calculation for different phase angles, then one single calculator node is sufficient, because the same can be reused for all angles and recalculated.

2) Next we have to calculate the particle velocity field, because it is needed for the Intensity calculations. The equations we use for this are:
Notation:
Eq_2.png
Eq_2.png (95.96 KiB) Viewed 8500 times
The velocity can be calculated from the pressure in Paraview by applying the filter “Gradient of Unstructured Dataset” to the input pressure dataset. Unfortunately this filter is not very professional, because the gradient values at the end boundaries of the tube are about -313 instead of 628, which is half of the correct value. This bug has been mentioned by mzenker in another thread where he also wrote that he has informed the Paraview mailing list about it. Any updates about this? My version of Paraview by the way is 4.4.

Only 3 attachments are allowed/post, therefore this text will continue below >

Zoltan
Last edited by Z.Losonc on 08 Jan 2018, 16:29, edited 5 times in total.
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by Z.Losonc »

...continued from above:

Here is the gradient calculated by Paraview:
tube_Paraview_Gradient_p2.png
tube_Paraview_Gradient_p2.png (73.73 KiB) Viewed 8695 times
This forces us to generate the velocity field using ElmerSolver, which is actually the best way to do it even if the Paraview filter gets fixed. The best would be to calculate the velocity as flux directly using the “FluxSolver”. In this case we would have to generate two new datasets, one for the real part and one for the imaginary part of the velocity. In this case the
Flux Coefficient = 1/(Density* Angular Frequency) which is 3.85e-4 here.
Re(v)=v1=-grad(Pressure Wave 2)/(Density*Angular Frequency)
Im(v)=v2=grad(Pressure Wave 1)/(Density*Angular Frequency)

Unfortunately I was not able to get Elmer 7 (on Windows) to produce the expected flux.
Does anybody know how to make this work?

However, the calculation of gradient did work and the results are correct at the boundaries as well. This way Elmer creates 4 datasets:
pressure wave 1
pressure wave 2
pressure wave 1 grad
pressure wave 2 grad
and the velocity has to be calculated in Paraview using these calculator formulas:
v1=-pressure wave 2 grad*3.85e-4
v2=pressure wave 1 grad*3.85e-4
The real part of velocity field v1 (x component) looks like this:
tube_Re_v.png
tube_Re_v.png (58.25 KiB) Viewed 8695 times
Here is the .sif file:

Code: Select all

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "."
  Include Path ""
  Results Directory ""
End

Simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Steady state
  Steady State Max Iterations = 10
  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 = Flux and Gradient 1
  Target Variable = Pressure Wave 1
  Procedure = "FluxSolver" "FluxSolver"
  Calculate Grad = True
  !Calculate Grad Magnitude = True
  !Discontinuous Galerkin = Logical True
  Exec Solver = After Simulation
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Linear System Solver = Direct
  Linear System Direct Method = Banded
End

Solver 2
  Equation = Flux and Gradient 2
  Target Variable = Pressure Wave 2
  Procedure = "FluxSolver" "FluxSolver"
  Calculate Grad = True
  !Calculate Grad Magnitude = True
  !Discontinuous Galerkin = Logical True
  Exec Solver = After Simulation
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Linear System Solver = Direct
  Linear System Direct Method = Banded
End

Solver 3
  Equation = Result Output
  Procedure = "ResultOutputSolve" "ResultOutputSolver"
  Output File Name = tube
  Output Format = Vtu
  Exec Solver = After Simulation
End

Solver 4
  Equation = Helmholtz Equation
  Procedure = "HelmholtzSolve" "HelmholtzSolver"
  Variable = -dofs 2 Pressure Wave
  Exec Solver = Before Simulation
  Stabilize = True
  Bubbles = False
  Lumped Mass Matrix = False
  Optimize Bandwidth = True
  Steady State Convergence Tolerance = 1.0e-5
  Linear System Solver = Direct
  Linear System Direct Method = Banded
End

Equation 1
  Name = "Helmholtz"
  Angular Frequency = 2155
  Active Solvers(4) = 1 2 3 4
End

Material 1
  Name = "Air (room temperature)"
  Heat Conductivity = 0.0257
  Heat Capacity = 1005.0
  Density = 1.205
  Relative Permittivity = 1.00059
  Viscosity = 1.983e-5
  Sound speed = 343.0
  Heat expansion Coefficient = 3.43e-3
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

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
When using an iterative solver as an alternative in the FluxSolver, which is based on a sample from mzenker the solution was OK, but the solver log contained messages that it did not converge.
Does anybody know why did this happen, and how to make it converge with iterative solver?

Code: Select all

Solver 1
  Equation = Flux and Gradient 1
  Target Variable = Pressure Wave 1
  Procedure = "FluxSolver" "FluxSolver"
  Calculate Grad = True
  !Discontinuous Galerkin = Logical True
  !Calculate Grad Magnitude = True
  Calculate Flux = True
  Calculate Flux Magnitude = True
  Flux Variable = String "Velocity 1"
  Flux Coefficient = Heat Conductivity
  Exec Solver = After Simulation
  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
End
Another problem I have encountered is that this project which was created in Elmer 7 (Windows) did not work in Elmer 6.2 (CAELinux 2013), and vice versa. Is there compatibility between versions? How to convert projects from one version to another?

After calculating the velocity datasets in Paraview one can create a velocity field animation the same way we did for the pressure using this equation (here fi=90deg):

v(fi)=v1*cos(90*3.14/180)-v2*sin(90*3.14/180)

The dataset for the absolute pressure can be calculated from the original Elmer dataset as:

p_abs=sqrt(pressure wave 1^2+pressure wave 2^2)
tube_p_abs.png
tube_p_abs.png (22.33 KiB) Viewed 8695 times
The radiation intensity can be calculated from the original pressure datasets plus the derived velocity datasets as:

Re_I=0.5*(pressure wave 1*v1+pressure wave 2*v2)
Im_I=0.5*(pressure wave 2*v1-pressure wave 1*v2)

Only 3 attachments are allowed/post, therefore this text will continue below >

Zoltan
Last edited by Z.Losonc on 27 Dec 2017, 01:02, edited 2 times in total.
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by Z.Losonc »

...continued from above:

Here is an image of the real radiation intensity field (x component):
tube_Re_Irms.png
tube_Re_Irms.png (22.55 KiB) Viewed 8695 times
The sound power that passes through the input boundary can be calculated in two ways:
a) From the pipeline node Re(I) create a plane slice at
origin: 0.000001, 0, 0 (when x=0 selected, got only half the input boundary; thus 0.000001)
normal: 1,0,0
From this node apply the filter “Surface Flow” > Select Input Vectors = Re_I
This will calculate the Pin=0.0943799 W which can be viewed in a SpreadSheetView Window.
tube_ParaView_Pipeline.png
tube_ParaView_Pipeline.png (20.62 KiB) Viewed 8695 times
tube_Pin_SurfaceFlow.png
tube_Pin_SurfaceFlow.png (9.09 KiB) Viewed 8695 times
b) The other (more complicated) way to get the same value is to apply the “Generate Surface Normals” filter on the “Input surface” node. Then use a calculator on this node with the formula: Re_I.Normals which will calculate the dotproduct of the Intensity vector and the surface normal vector. Finally apply an “Integrate Variables” filter on this last node, and you will get the same Pin value as above.

Finally the sound pressure level is calculated from the original pressure field as:

SPL=10*log10(0.5*(pressure wave 1^2+pressure wave 2^2)/2.0E-5^2)

Only 3 attachments are allowed/post, therefore this text will continue below >

Zoltan
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by Z.Losonc »

...continued from above:

And here is how it looks like:
tube_SPL.png
tube_SPL.png (23.56 KiB) Viewed 8695 times
I am looking forward to hear your feedback about this analysis, whether it is correct or there are any mistakes to be corrected. It would be great if somebody could explain (and show) how to generate at least the velocity field (but preferably also the p_abs, SPL, and radiation intensity fields) using ElmerSolver.

The project folder is attached below, but due to file size limitation I could not include all the necessary files. Therefore, after unzipping the project folder move the Simple_Pipe_1.mphtxt mesh file into it, and run the solver. If you are a beginner and want to use the ElmerGUI (like I do), then open the project (this will automatically overwrite the tube.sif with a basic one that doesn’t generate gradients). Now open the manually written backup tube.sif0 in a text editor > copy all the text > open tube.sif in another text editor > delte all the text > paste the copied text > save tube.sif.

Next, (without generating Sif !) run the solver. You can modify the sif either in ElmerGUI > Sif > Edit, or in an external text editor; save, and rerun the solver to experiment with different settings. But if you close ElmerGUI and restart it for whatever reason, then you will have to redo the copy/paste operation from the tube.sif0 backup as described above. If you don’t want to loose your modified sif, then make a backup copy that serves the function of tube.sif0.

The file tube.pvsm is the Paraview state file that you can load to see the ready made postprocessing.
Only 3 attachments are allowed/post, therefore the rest of the files will be attached below >

Zoltan
Attachments
Simple_Pipe_1.mphtxt
(815.63 KiB) Downloaded 405 times
Acoustic_Tube_Postprocessing.7z
(157.71 KiB) Downloaded 414 times
Last edited by Z.Losonc on 22 Dec 2017, 12:28, edited 1 time in total.
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by Z.Losonc »

...continued from above:

And finally two more attachments below.

Cheers,
Zoltan
Attachments
tube.7z
(14.85 KiB) Downloaded 424 times
tube0001.7z
(362.64 KiB) Downloaded 412 times
mzenker
Posts: 1999
Joined: 07 Dec 2009, 11:49
Location: Germany

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by mzenker »

Hi,

I am no acoustics expert, so I cannot say if your results are correct. But I would recommend to update your Elmer under Windows and Linux to the actual version 8. Elmer 7 and even more Elmer 6 are outdated, and I would not invest time in trying to solve problems you get with old software versions.

Matthias
annier
Posts: 1168
Joined: 27 Aug 2013, 13:51
Antispam: Yes

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by annier »

Hi Zoltan,
1. I second Matthias in recommending you to use the latest Elmer versions . Peter, Juha and other members of Elmer team can answer your queries meticulously provided you are using latest versions of Elmer.

2. If you are using Ubuntu 16.04 , installation of Elmer by launchpad is the easiest and complete way of utilizing the latest features of the multiphysics software. Peter loves several queries on Elmer Solver and answers them frequently.

3. Nice presentation.


4. The best way to learn Elmer is to look at the source codes of the Elmer Solver as provided in github repository and all the tutorials/examples direct to the source codes.


Yours Sincerely,
Anil Kunwar
Anil Kunwar
Faculty of Mechanical Engineering, Silesian University of Technology, Gliwice
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by Z.Losonc »

Hi Matthias and Anil,

Thank you for your comments and suggestions.

I am working on the compilation of the latest version of Elmer on Linux Mint 17, which will take some time, because I will have to reorganize my hard drive first to make space for an extra partition.

The reason for me using the older 6.2 and 7 versions is that they were installed some years ago out of curiosity, just to see if I could get it work for acoustics. There were no ready made tutorials in acoustics at the time that I could simply load, run and see useful result; and did not have the time and patience to dig deeper into the documentation to figure things out. Then this idea faded into the background and got forgotten. But now Joe’s ready made project got me started on this again.

In mean time I have found the problem, and successfully computed the flux (producing the acoustic velocity in this case) with the FluxSolver on Elmer 7, so it looks like this feature is not version dependent. I will explain this in detail in my next post after re-running the project on the latest Elmer version as you have suggested.

Regards,
Zoltan
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by Z.Losonc »

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
Z.Losonc
Posts: 40
Joined: 21 Dec 2017, 16:05
Antispam: Yes

Re: Acoustics – Basic Postprocessing With ElmerSolver & Paraview

Post by Z.Losonc »

An update: found out how to define an arbitrary Flux Coefficient for the FluxSolver and assign a value to it via a formula.

In the “Material 1” section define the arbitrary variable (or constant) Cv either by assigning a known number to it like:

Cv = Real 3.85e-4

or by assigning a value to it by a MATC expression, which should use known parameters as arguments defined elsewhere:

Cv = Real MATC "1/(Density * Angular Frequency)"

Finally, replace the default form of Flux Coefficient declaration in the FluxSolver sections (Flux Coefficient = String "Heat Conductivity") with:

Flux Coefficient = String "Cv"

(Cv is an arbitrarily chosen name for the Flux Coefficient parameter; replace it with your own chosen name)

Zoltan
Post Reply