microstructure_history_MFIX.m

%% ------------------------------------------- Script to gather relevant simulation data from MFIX history ------------------------------------------- %%

% Rajath Kantharaj
% Purdue University
% CTRC TIMs project, January 2019 - December 2020

% -------------------- Inputs -------------------- %
% 1. FIND_flag = 'BLT' or 'Pressure': Capture end state of simulation
% based on BLT or squeeze pressure provided by user
% 2. blt_find -- Capture state of simulation with desired BLT
% 3. pressure_find -- Capture state of simulation with desired pressure
% 4. Path_File -- path to (multiple) folders containing 'wall_dynamics.txt'
%    file. Specify as a cell vector
% 5. DT_VTK_vec -- time interval for writing VTK files in the simulation
% 6. DTSOLIDS -- DEM time step in the simulation (usually 1e-3 or 0.5e-3)

clc; clearvars; close all
%% ------------------------------------------- INPUTS ------------------------------------------- %%
%% Match BLT or pressure?
FIND_flag = 'Pressure'; % 'BLT'
blt_find = 2.9687e-04; % XRCT squeezed TIM BLT
pressure_find = 150e3; % units: Pa (for reference,  25 PSI ~ 172.4 kPa)

%% Path, file names and other inputs
Path_File = cellstr(input('Enter the path to MFIX output files \n','s'));

%% Define time steps and other parameters
DT_VTK_vec = input('Enter the time step for writing output VTK files \n'); % 0.001; %[0.001, 0.001, 0.00025]; % Write time interval of VTK files --> see MFIX input file
DTSOLIDS = input('Enter the DEM simulation time step \n') ; %1.0408e-9; % ; % seconds - DEM solids time step --> see output from MFIX 

% ------------------------------------------------------------------------------------------------ %

%% ------------------------------------------- Do not modify code below this line ------------------------------------------- %%

write_freq = 50; % value of "calls_write" variable from MFIX in "usr1_des.f" module
WALL_file = 'wall_dynamics.txt'; % file containing pressure, force and other crucial information

%BLT_file = 'blt_plate_pos.txt';
%nDp = 3; % Find microstructure info for BLT = nDp * max(Dp);

% Loop through and compute
for ii=1:numel(Path_File)
    path_file = Path_File{ii};
    DT_VTK = DT_VTK_vec(ii);
        
    % ---- Step 1 ----
    % Plot BLT from wall_dynamics TXT file & check time stamp for BLT of
    % interest, "blt_find"
    B = importdata([path_file '\' WALL_file]);
    blt_sim = B.data(:,2); % BLT from simulation
    prss = B.data(:,4); % predicted pressure
    
    tol_blt = 1e-3; % Does the tolerance need to be this tight?
    
    if strcmpi('BLT',FIND_flag)
        indx = find((blt_sim - blt_find)./blt_sim < tol_blt, 1);
        fprintf('MFIX: Bond line thickness BLT = %.3e m \n', blt_sim(indx));
    elseif strcmpi('Pressure', FIND_flag)
        indx = find(abs(prss - pressure_find)./prss < tol_blt, 1);
        fprintf('MFIX: Squeeze pressure is P = %.2f kPa \n', prss(indx)/1e3)
        
        fprintf('BLT = %3.2f um \n', blt_sim(indx)*1e6)
    end
    
    % ---- Step 2 ----
    indc = [2:size(blt_sim,1)]' - 1;
    time_wall = [0; write_freq*DTSOLIDS*indc];
    real_time = time_wall(indx);
    
    time_stamp_VTK = floor(real_time/DT_VTK);
    
    cp = strsplit(path_file,'\'); 
    cp = cell2mat(cp(end));
    fprintf('The vtp file number for file %s is %d\n', cp, time_stamp_VTK)

end