English
Log in
English

How to batch process your OnScale Sweep Results

By Cyprien Rusu 16 September 2020

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.

PZT PillarsOur 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:

Sweep

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

Piezoelectric

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

Design Variables

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:

Parametric Sweep

After clicking on run, you will see the 400 simulations starting to run in parallel on the cloudSweep Results

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

OnScale Analyst

Then save it under the “revinp” format. We can call it “process.revinp” for example.

sweep result

Make sure that you save that script in the folder where results are located:

review script

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”

sweep data manager

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:

script

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:

sweep process

We will now get the following kpi_report.csv with the following content:

KPI report

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:

Post Processor

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:

Sweep Results

If we uncomment the wndo command, we observe that the time curve now contains lesser harmonics:

Sweep Results

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):
sweep results

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

graph command

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

3D composite

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:

KPI

Run Review:

And you will get all the data for the 400 simulations post-processed into your excel file!

simulation results

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

MATLAB results

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.

Cyprien Rusu
Cyprien Rusu

Cyprien Rusu is our Director of Engineering for Asia at OnScale. He has a extensive background in FEM, Technical Marketing, Sales and Support. Cyprien recieved his MS in Civil Engineering from Tsinghua University. At OnScale he is a trusted advisor for our client base in Asia, while creating and providing OnScale training.