One of the great features of OnScale is that it is easy to launch a large parametric sweep on the cloud and get thousands of results at the same time.
That’s impressive!
In this previous blog post, we discussed how to do that and why: Simulation and optimization of piezocomposite transducers with a practical example
But the great question we are going to tackle here is… how do you process those thousands of data files generated?
You can’t do that by hand that’s for sure!
That’s when the power of review post-processing automation comes into play.
In this article, I will show you step by step how to create an automation script that will do the following:
-
Open and read thousands of result files
-
Process and transform data stored within those files
-
Save a summary of what interests us into a simple csv file
We will do all that by building step by step a review file which will perform all the important functions.
1- Presentation of the model
Let’s look now at the following specific problem involving a 1 MHz piezocomposite sensor operating into water and being used in pulse-echo mode against a steel target.
Our objective of this analysis is to obtain a high sensitivity and bandwidth when operating in pulse echo mode.
For this exercise, we will make some initial design decisions:
-
We are going to use a soft PZT with a high coupling coefficient, because this will provide us a high sensitivity and bandwidth
-
We are also going to use a hard epoxy filler
This is what we take for given and we won’t change
Here’s a summary of the materials chosen:
Detailed material properties are provided in the model file “composite.prjmat”. Those properties are also available in OnScale Material database
2- Running the Simulation Parametric Sweep
We want to maximize sensitivity and bandwidth and we want to keep the center frequency around 1 MHz
The parameters that we are going to vary are:
-
The matching layer thickness (between 0.45 and 0.75 mm)
-
And we will vary the composite volume fraction between 40-60%
(This is a typical range for high bandwidth pulse-echo mode transducers)
We will monitor the 3 following KPIs:
-
the sensitivity
-
the bandwidth
-
the center frequency
We will now use OnScale to perform a parametric sweep of those 2 design variables.
In total, we will get 20 values of matching layer * 20 values of Composite volume fraction, for a total of 20*20 = 400 simulations!
Let’s launch the simulation in OnScale and choose the following definition for the 2 parametric design variables:
After clicking on run, you will see the 400 simulations starting to run in parallel on the cloud
Wait until the end.
The calculation time of those 400 simulations takes around 7 min in total and around 100 CH !
When this is over, download the results from the cloud.
Results we are interested in are included in a lightweight composite_disk.flxhst Time history file, so the download process is also fast.
3- Batch Post-processing the results with Review
Now that we got all the results download in our local computer, let’s see how to post-process those results using a review script.
Review is using the Symbol language as a script language. If you want a detailed introduction to Symbol, you can watch this video.
Let’s start by opening a new file in OnScale Analyst Mode and
Then save it under the “revinp” format. We can call it “process.revinp” for example.
Make sure that you save that script in the folder where results are located:
Step 1 – Setting up the post-processing loop
Paste the following code in it.
process.revinp
c Multithreading mp omp * * c Define variables symb nfiles = 400 text tfile = 'composite_disk.flxhst' text tdir = 'composite_disk-20200324-165920' c Loop through files do loopi I 1 $nfiles 1 c Add the code to process the files here... end$ loopi stop
Explanation of this code:
1- The “mp omp * *” command is added so review will use multithreading to process results (this will be much faster!)
2- Let’s define a few variables:
-
the nfiles variables represents the number of files that we will post-process
-
the variables tfile and tdir are text strings that will be used to read the results
3- Let’s now define a loop which will run from 1 to $nfiles times (400 times in total):
-
We open the loop using “do loopi I 1 $nfiles 1”
-
We close the loop using “loopi”
Let’s now see how to read content from files and save it to one unique csv file
Step 2 – Reading the result files and writing to one unique csv result file
Now, let’s process the data included within the result files!
Let’s add the following code within the loop.
...
c Report progress
symb #msg 1
Reading file $I of $nfiles ...
c Read data
text tpath = '$tdir/$I/$tfile'
read f1 $tpath
symb #get { valmax } curvmax f1 1
symb #msg 1
$valmax
...
This code will help us to test if we are able to read the content of the result files and print some data in the onscale console.
Explanation of the code:
1- The command “symb #msg 1” is the command used to print a message of 1 line to the console
2- The following line “Reading file $I of $nfiles …” will then be printed to the console and the 2 variables $I and $nfiles will be automatically replaced
3- The string variable tpath will contain the path to the result file.
note the peculiar syntax:
text tpath = '$tdir/$I/$tfile'
This will become for example the following path when the result number 1 of the loop is processed!
text tpath = 'composite_disk-20200324-165920/1/composite_disk.flxhst'
4- Now, let’s read the data and save it in a file called f1 with the command “read f1 $tpath“
5- We now have access to the data within the result file which have been saved into the file f1
If we stopped the algorithm with a term command at every iteration, we would be able to see the file f1 appearing in the Data Manager. The time curve data we are interested in are saved in the array “time”
6- Let’s save the maximum of the time curve stored in first position in each f1 file in another variable valmax
for that we use the following function curvmax
symb #get { valmax } curvmax f1 1
This code has to be understood like that:
Create a new valmax variable and assign to it the maximum value of the first curve stored in f1.
Note that valmax is a real, so the first letter of the variable should be in the range A-H and O-Z (Fortran convention). All the names starting by a letter within I-N are considered as integers!
7- Let’s now print again this value to the console with:
symb #msg 1 $valmax
If you now run review for 10 iterations, you should see the following displayed in the console:
This means that the script works!
Now, instead of the console, let’s add some additional code to print out the result to a file
...
c Write report file header
symb #msg 1 > 'kpi_report.csv'
Model, valmax
c Loop through files
do loopi I 1 $nfiles 1
c Report progress
symb #msg 1
Reading file $I of $nfiles ...
c Read data
text tpath = '$tdir/$I/$tfile'
read f1 $tpath
symb #get { valmax } curvmax f1 1
c Write to screen and file
symb #msg 1
$I, $valmax
symb #msg 1 >> 'kpi_report.csv'
$I, $valmax
end$ loopi
...
Explanation:
1- Right before the opening of the loop we added a file header that will create a csv file:
c Write report file header symb #msg 1 > 'kpi_report.csv' Model, valmax
This is exactly like printing to the console, except that we are now adding a > ‘kpi_report.csv’ which will redirect the content to a file. The “>” means that we want to generate a new file. Thus every file with the name ‘kpi_report.csv’ which was present in the folder before will be deleted and written over by this new message print.
2- Inside the loop, but right before the closing end$, we added the following code
symb #msg 1 >> 'kpi_report.csv' $I, $valmax
This means that the message $I, $valmax will be appended to the file kpi_report.csv at every iteration.
Note that now we use “>>” to signify that we want to append to the existing file!
If you run review again for 10 iterations:
We will now get the following kpi_report.csv with the following content:
Note: A csv file is a simple text file. You can also open it with a simple text editor!
Step 3 – Obtaining meaningful KPIs
Now, let’s go deeper and let’s actually calculate some useful KPIs using Review!
Let’s open first one of the composite_unit_cell.flxhst result files from within OnScale Post-process GUI:
We see that we have 5 time history curves stored in it.
What interests us is the 3rd curve which contains the full time signal.
We can then use the Plot Control> Markers to spot what is the time range that really interest us.
Ok, let’s get back to the script now!
The following script placed the read command will help us to start processing the data within the file f1.
symb nt = 1500
symb tb = 1.4e-5
symb ts = 5.e-9
...
read f1 $tpath
c Get pulse
make
file f1w cler
wndo hann righ 0. 0.e-6 20.e-6
time $nt $tb $ts
curv { f1 3 }
end
grph
plot f1w 1
term
...
The make command is used here to create a new file called f1w which will contain an extracted part of the time curve.
The grph command is used to plot the content of the new generated curve within the file f1w
The term command allows us to pause the run after each loop iteration to leave us some time to review the graphics.
the file command opens the file. Note the “cler” option which is important here. It is used to clear all records previously stored in local file f1w.
the wndo command is used to apply a window function which will limit the frequency harmonics included within the final result. It is used here to remove lower harmonics susceptible to cause problem and unwanted noise during a conversion to the frequency domain.
the time command will limit the part of the time domain processed. Note that we added 3 new variables $nt $tb $ts which represent respectively the number of time steps, the beginning time and the size of the time increment.
the curv command is here to provide the actual data that will be processed in input. In this case, the 3rd time history curve.
If we run review now, we should get the following:
If we uncomment the wndo command, we observe that the time curve now contains lesser harmonics:
Let’s now change that code to calculate the FFT of the signal:
symb nt = 1500
symb tb = 1.4e-5
symb ts = 5.e-9
symb fmin = 0.e6
symb fmax = 2.e6
symb nf = 5000
symb finc = ( $fmax - $fmin ) / $nf
...
c Get pulse
make
file f1w cler
wndo hann righ 0. 0.e-6 20.e-6
time $nt $tb $ts
curv { f1 3 }
end
c Calculate fft
freq
file f2 cler
wndo hann righ 1
type amp
fft f1w 1
end
c Group fft data
make
file f3
freq $nf $fmin $finc
curv { f2 1 }
end
grph
plot f3 1
term
Here, the important command is the command freq which handles the conversion of the time data into frequency data:
c Calculate fft freq file f2 cler wndo hann righ 1 type amp fft f1w 1 end
The sub-command type amp is used to specify we want data saved in amplitude/phase format rather than real/imaginary.
The sub-command freq fft fw1 1 specifies that the curve to transform with FFT is the curve stored in the 1st place of the f1w file.
We then post-process those data a second time in order to limit the frequency range.
symb fmin = 0.e6
symb fmax = 2.e6
symb nf = 5000
symb finc = ( $fmax - $fmin ) / $nf
...
c Group fft data
make
file f3
freq $nf $fmin $finc
curv { f2 1 }
end
Note the 4 new variables added at the top of the review file.
Ok, Now just remain to plot the frequency domain curve:
grph plot f3 1
Run Review and this is what you should be getting (at each iteration):
We can see that the resolution of that curve isn’t that great!
That’s because we have limited the number of time samples to 1500. Let’s change nt to 5000 time samples.
symb nt = 5000
Let’s now modify again the script and replace the current grph command by the following
symb #get { rdelx } curvdelx f1w 1
symb #get { nval } curvnval f1w 1
symb #msg 1
Nsamples = $nval, Ts = $rdelx
grph
nvew 2 1
pset wndo off
plot f1w 1
pset wndo $fmin $fmax
plot f2 1
This will upgrade the plot to visualize the time and the converted frequency signals at the same time.
-
the function curvdelx returns the difference between steps on the x axis
-
the function curvnval returns the number of x-y values in the local file
Ok, now we just have to calculate the important KPIs we are interested in.
Let’s add the following script right before the end of the loop
c Group time data
make
file f4
curv { f1w 1 }
end
c Get centre frequency and amplitude
symb #get { fc amp } curvmax f3 $I $fmin $fmax
c Get BW
symb amp2 = $amp / 2.
symb #get { XL } curvcros f3 $I $amp2 GL 1 $fmin $fc
symb #get { XR } curvcros f3 $I $amp2 GR 1 $fc $fmax
symb bandwidth = ( $GR1 - $GL1 ) / 1.e6
symb fcmhz = ( $GL1 + $bandwidth / 2. ) / 1.e6
symb FBW = 100. * $bandwidth / $fcmhz
c Write to screen and file
symb #msg 1
$I, $fcmhz, $amp, $FBW
symb #msg 1 >> 'kpi_report.csv'
$I, $fcmhz, $amp, $FBW
1- the function curvmax extracts the center frequency fc and the corresponding amplitude from the FFTed frequency domain curve stored in f3.
2- We then calculate the Bandwidth
3- And finally we calculate the KPIs we want:
Run Review:
And you will get all the data for the 400 simulations post-processed into your excel file!
Step 3 – Saving all data to CSV Files
We will now export the data stored into the f3 and f4 files using the data cddo command to create csv files:
c Write out data data cddo f3/data fft.csv cddo f3/time freq.csv cddo f4/data amp.csv cddo f4/time time.csv end
If you open the folder where your model is saved, you should now find the 5 following Excel CSV files:
Those files contain all Time-Amplitude and Frequency-FFT data for all the 400 simulations!
We will use those files as an input to a Matlab script which will allow us to post-process the results graphically.
Here is the final script used:
c Multithreading
mp omp * *
c Define variables
symb nt = 5000
symb tb = 1.4e-5
symb ts = 5.e-9
symb fmin = 0.e6
symb fmax = 2.e6
symb nf = 1000
c Calcs
symb finc = ( $fmax - $fmin ) / $nf
symb nfiles = 400
text tfile = 'composite_disk.flxhst'
text tdir = 'composite_disk-20200324-165920'
c Write report file header
symb #msg 1 > 'kpi_report.csv'
Model, Fc, Amp, FBW
c Loop through files
do loopi I 1 $nfiles 1
c Report progress
symb #msg 1
Reading file $I of $nfiles ...
c Read data
text tpath = '$tdir/$I/$tfile'
read f1 $tpath
c Get pulse
make
file f1w cler
wndo hann righ 0. 0.e-6 20.e-6
time $nt $tb $ts
curv { f1 3 }
end
c Calculate fft
freq
file f2 cler
wndo hann righ 1
/* time pad $npad
type amp
fft f1w 1
end
c Group fft data
make
file f3
freq $nf $fmin $finc
curv { f2 1 }
end
symb #get { rdelx } curvdelx f1w 1
symb #get { nval } curvnval f1w 1
symb #msg 1
Nsamples = $nval, Ts = $rdelx
grph
nvew 2 1
pset wndo off
plot f1w 1
pset wndo $fmin $fmax
plot f2 1
c term
c Group time data
make
file f4
curv { f1w 1 }
end
c Get centre frequency and
c symb nrec = $I * 2 - 1
symb #get { fc amp } curvmax f3 $I $fmin $fmax
c symb fcmhz = $fc / 1.e6
c Get BW
symb amp2 = $amp / 2.
symb #get { XL } curvcros f3 $I $amp2 GL 1 $fmin $fc
symb #get { XR } curvcros f3 $I $amp2 GR 1 $fc $fmax
symb bandwidth = ( $GR1 - $GL1 ) / 1.e6
symb fcmhz = ( $GL1 + $bandwidth / 2. ) / 1.e6
symb FBW = 100. * $bandwidth / $fcmhz
c Write to screen and file
symb #msg 1
$I, $fcmhz, $amp, $FBW
symb #msg 1 >> 'kpi_report.csv'
$I, $fcmhz, $amp, $FBW
end$ loopi
c Write out data
data
cddo f3/data fft.csv
cddo f3/time freq.csv
cddo f4/data amp.csv
cddo f4/time time.csv
end
stop
4- Graphical Visualization of the results in Matlab
In order to visualize the results in Matlab, we will use the following Matlab script.
batch_plot_script_link.m
clear;
close all;
clc;
% Inputs
V1Name = 'Matching Thickness (mm)';
V1Beg = 0.45;
V1End = 0.75;
% V1Beg = 0.52;
% V1End = 0.58;
NV1 = 20;
V2Name = 'Volume Fraction (%)';
% V2Beg = 0.45;
% V2End = 0.55;
V2Beg = 40;
V2End = 60;
NV2 = 20;
FMin = 0;
FMax = 2;
TMin = 10;
TMax = 30;
% Create vectors
V1Vect = linspace(V1Beg,V1End,NV1)';
V2Vect = linspace(V2Beg,V2End,NV2)';
[V1Grid,V2Grid] = meshgrid(V1Vect,V2Vect);
PlotSub = [1,1];
PlotInd = sub2ind([NV2,NV1],PlotSub(1),PlotSub(2));
% Pull results
Freq = csvread('freq.csv')';
Freq = Freq(1:end-1);
Nf = length(Freq);
Time = csvread('time.csv')';
Time = Time(1:end-1);
Nt = length(Time);
FFT = csvread('fft.csv')';
FFT = FFT(1:end-1,:)*1e6;
Volt = csvread('amp.csv')';
Volt = Volt(1:end-1,:)*1e3;
KPI = csvread('kpi_report.csv',1);
% Create KPI Matrices
Fc = reshape(KPI(:,2),NV2,NV1);
Amp = reshape(KPI(:,3),NV2,NV1)*1e6;
FBW = reshape(KPI(:,4),NV2,NV1);
JobNum = reshape(KPI(:,1),NV2,NV1);
% Get max amp for plotting
AmpMax1 = 1.05 * max(max(FFT(1:end-1,:)));
AmpLim1 = [0 AmpMax1];
AmpMin2 = 1.05 * min(min(Volt(1:end-1,:)));
AmpMax2 = 1.05 * max(max(Volt(1:end-1,:)));
AmpLim2 = [AmpMin2 AmpMax2];
% Create figure and setup datatip
f = figure('Name','OnScale Post Processing','WindowStyle','normal','NumberTitle','off',...
'units','normalized','outerposition',[0.05 0.075 0.9 0.9]);
cursorMode = datacursormode(f);
set(cursorMode,'enable','on','UpdateFcn',@gettippos);
% Colour codes
CBlue = [0 0.4470 0.7410];
CGreen = [0.4660 0.6740 0.1880];
% Plot fft
subplot(2,2,1);
% figure;
% yyaxis left;
% hPlot = plot(Freq/1e3,FFT(:,PlotInd));
hline1 = line(Freq/1e6,FFT(:,PlotInd));
xlim([FMin FMax]);
xlabel('Frequency (MHz)');
ylabel('FFT Magnitude');
ylim(AmpLim1);
hax1a = gca;
hax1a.XColor = CBlue;
hax1a.YColor = CBlue;
hax1a_pos = hax1a.Position;
hax1b = axes('Position',hax1a_pos,...
'XAxisLocation','top',...
'YAxisLocation','right',...
'Color','none');
hline2 = line(Time*1e6,Volt(:,PlotInd),'Parent',hax1b,'Color',CGreen);
hax1b.XColor = CGreen;
hax1b.YColor = CGreen;
xlabel('Time (us)');
ylabel('Voltage (mV)');
xlim([TMin TMax]);
ylim(AmpLim2);
% Plot sensitivity & add datatip
hax1 = subplot(2,2,2);
hSurf = surf(V1Vect,V2Vect,Amp,'FaceColor','interp','FaceLighting','gouraud');
QHand = gca;
set(QHand,'Xdir','reverse','Ydir','reverse');
hDatatip = cursorMode.createDatatip(hSurf);
set(hDatatip,'Position',[V1Beg,V2Beg]);
title('Sensitivity');
xlabel(V1Name);
ylabel(V2Name);
zlabel('Sensitivity (mV)');
xlim([V1Beg V1End]);
ylim([V2Beg V2End]);
% Plot frequency
hax2 = subplot(2,2,3);
surf(V1Vect,V2Vect,Fc,'FaceColor','interp');
FsHand = gca;
set(FsHand,'Xdir','reverse','Ydir','reverse');
title('Frequency');
xlabel(V1Name);
ylabel(V2Name);
zlabel('Frequency (Hz)');
xlim([V1Beg V1End]);
ylim([V2Beg V2End]);
% Plot bandwidth
hax3 = subplot(2,2,4);
surf(V1Vect,V2Vect,FBW,'FaceColor','interp');
FpHand = gca;
set(FpHand,'Xdir','reverse','Ydir','reverse');
title('Fractional Bandwidth');
xlabel(V1Name);
ylabel(V2Name);
zlabel('Fractional Bandwidth (%)');
xlim([V1Beg V1End]);
ylim([V2Beg V2End]);
% caxis([0 150]);
% Link rotation on plots
% hlink = linkprop([hax1,hax2,hax3],{'CameraPosition','CameraUpVector'});
hlink = linkprop([hax1,hax2,hax3],{'View'});
% rotate3d on;
% Copy data to datatip to be accessible from update function
hSurf.UserData.FFT = FFT;
hSurf.UserData.Volt = Volt;
hSurf.UserData.V1Vect = V1Vect;
hSurf.UserData.V2Vect = V2Vect;
hSurf.UserData.V1Name = V1Name;
hSurf.UserData.V2Name = V2Name;
hSurf.UserData.hline1 = hline1;
hSurf.UserData.hline2 = hline2;
% hSurf.UserData.PlotYLim1 = AmpLim1;
% hSurf.UserData.PlotYLim2 = AmpLim2;
% Datatip update function
function TipText = gettippos(~,hTarg)
% Get position
TipPos = get(hTarg,'Position');
% Get indices
NV1 = length(hTarg.Target.UserData.V1Vect);
NV2 = length(hTarg.Target.UserData.V2Vect);
PlotSub(1) = find(hTarg.Target.UserData.V2Vect == TipPos(2),1,'first');
PlotSub(2) = find(hTarg.Target.UserData.V1Vect == TipPos(1),1,'first');
PlotInd = sub2ind([NV2,NV1],PlotSub(1),PlotSub(2));
% Update plot
set(hTarg.Target.UserData.hline1,'YData',hTarg.Target.UserData.FFT(:,PlotInd)');
set(hTarg.Target.UserData.hline2,'YData',hTarg.Target.UserData.Volt(:,PlotInd)');
% set(hTarg.Target.UserData.hPlot,'YLim',hTarg.Target.UserData.PlotYLim);
% Create tip text
TipText = {[hTarg.Target.UserData.V1Name ': ' num2str(TipPos(1))]; ...
[hTarg.Target.UserData.V2Name ': ' num2str(TipPos(2))]; ...
['Sensitivity: ' num2str(hTarg.Target.ZData(PlotSub(1),PlotSub(2)),3)]};
end
The script has to be placed in the same folder than the results.
Then open this file with Matlab and click on Run
You then get a graphical representation in 3D of the KPIs.
The cursor can be moved for a specific KPI and the FFT and Voltage amplitude displayed then change dynamically.
… And this is how you batch process your sweep results in OnScale!
It may look complicated because this blog post details all the mechanic of the script so you can edit it or build your own, but in fact, you could also just reuse those scripts to generate the data and post-process with Matlab immediately even without going into the nuts and bolts of how it works.
Good luck and if you have any question, please leave them in the comments.
