Plotting residuals in openFOAM

There are many resources that will show you how to plot residuals while solving, for example look here. The core of most (or all..) of them is to filter the numeric values into gnuplot using few unix utilities like so:

set logscale y
set title "Residuals"
set ylabel "Residual"
set xlabel "Iteration"
plot "< cat log | grep 'Solving for Ux' | cut -d' ' -f9 | tr -d ','" title 'Ux' with lines,\
     "< cat log | grep 'Solving for Uy' | cut -d' ' -f9 | tr -d ','" title 'Uy' with lines,\
     "< cat log | grep 'Solving for Uz' | cut -d' ' -f9 | tr -d ','" title 'Uz' with lines,\
     "< cat log | grep 'Solving for omega' | cut -d' ' -f9 | tr -d ','" title 'omega' with lines,\
     "< cat log | grep 'Solving for k' | cut -d' ' -f9 | tr -d ','" title 'k' with lines,\
     "< cat log | grep 'Solving for p' | cut -d' ' -f9 | tr -d ','" title 'p' with lines
pause 1
reread

Notice that each line is piping the log file thru 4 different utilities. I think it is possible to make it simpler (but a bit slower.. if you don't care, continue reading), lets try it. The first thing to fly out is tr -d ',' which simply removes the semy-colon from the numbers. Turns out gnuplot can digest the semi-colons just fine. The next thing we can do is input the log file into grep instead of piping it, which should be faster, and we can look up the string 'for Ux' or whatever, dropping the 'Solving' part. Now we have (using only one of the lines as example)

grep 'for Ux' log | cut -d' ' -f9

That's already much shorter and faster. But, what we are doing here is looking for a fixed string in a file with many lines, taking all those lines and using just a number in column 9 of those lines. Turns out there is a utility that is optimized to do just that. It is called awk and like grep and cut it is installed on every unix/linux out there. That same line using awk now look like this:

awk '/for Ux/ {print $8}' log

Great, it is even shorter now and more readable. awk is looking for all lines with the string 'for Ux' and printing the 8th column. Apparently cut and awk are using different methods to calculate the columns, you can look it up if you like. Lets look at the gnuplot file now:

set logscale y
set title "Residuals"
set ylabel "Residual"
set xlabel "Iteration"
plot "< awk '/for Ux/ {print $8}' log" title 'Ux' with lines,\
     "< awk '/for Uy/ {print $8}' log" title 'Uy' with lines,\
     "< awk '/for Uz/ {print $8}' log" title 'Uz' with lines,\
     "< awk '/for omega/ {print $8}' log" title 'omega' with lines,\
     "< awk '/for k/ {print $8}' log" title 'k' with lines,\
     "< awk '/for p/ {print $8}' log" title 'p' with lines
pause 1
reread

Much better. Now when using a solver like simpleFoam, I usually solve for the pressure at least twice in every timestep. Doing so and plotting results in this:

We get twice (in this case) the number of pressure residuals as the rest of them, and it's alteranting between the first and second computation values, this is the reason for the thick band. Note that the computation is not converged here, it is just for demonstration.

OK, what is needed is another filter to take only every other pressure line. We can do it with awk as well, since it can do arithmatic operations and it knows the line numbers of the input file:

awk '/for Ux/ {print $8}' log | awk 'NR%2==0'

NR is the internal awk variable for line number, and we are checking for a reminder of zero when dividing by 2. If we need every third line we will divide by 3, etc. Why is there no {print} statement? Because the default action of awk is to print the whole line when the condition is true, and at this stage our whole line is just the number we are interested in.

Can we do better? Apparently we can... It turns out that gnuplot itself can filter blocks and values. We will use this to filter values (rows). It is done like so:

plot "< awk '/for p/ {print $8}' log" every 2::1 title 'p' with lines

What this means is 'plot every 2 values, starting at the 2nd value' (gnuplot numbering starts from 0). Look up gnuplot help for more info if interested. Now the same residuls plot looks like this:

Still not converged, but looking much better... This post is becoming long, so lets finish with one last improvement. Notice that the solver-log filename is appearing on each line of the plot command, but what if we are using another filename? Should we change all those lines? Of-coarse not, and here gnuplot user defined variables can help us. Consider this last and final listing:

log = "longlogfilename"
set logscale y
set title "Residuals"
set ylabel "Residual"
set xlabel "Iteration"
plot "< awk '/for Ux/ {print $8}' ".log title 'Ux' with lines,\
     "< awk '/for Uy/ {print $8}' ".log title 'Uy' with lines,\
     "< awk '/for Uz/ {print $8}' ".log title 'Uz' with lines,\
     "< awk '/for omega/ {print $8}' ".log title 'omega' with lines,\
     "< awk '/for k/ {print $8}' ".log title 'k' with lines,\
     "< awk '/for p/ {print $8}' ".log every 2::1 title 'p' with lines
pause 3
reread

log is now our variable which holds the filename, and we use a dot to concatenate it to the awk command.

If the pandemic continues much longer, I will probably find ways to make it shorter ;-)