Description of solvePbe

 solve a population balance (growth, nucleation, seeding, bulk states)

 Syntax: result = solvePopulationBalance(pbeDefinition)

 This routine essentially implements the method of characteristics for a
   set of input points. In this scope, the routine considers the
   limitation of the growth rate to the gemetrical boundaries (e.g.
   validity). It does also recognize dissolving particles.
 Not explicitely contained in the scheme of this function is the
   integration of the particle population to, e.g. the total particle
   volume. This functionality is implicitely provided by the definition of
   bulk properties in the PBE definition structure (pbeDef). However, this
   library follows the approach of Monte Carlo integral estimates which
   is, in particular, reflected by the function extractDataFromPbeResult.
 The support for multiple seeding and nucleation was cancelled because the
   overal scheme with Monte Carlo integral estimates does not yet cover
   these situations. It would, however, be straightforward to add this
   functionality. Your turn! ;)
 For input/output see the link below.

 Detailed documentation for inpu/output structures and option/value pairs
   for the function call are available, seperately, <a href="matlab: helpPbeSolver">here</a>.
   This documentation contains information on:
     <a href="matlab: help helpPbeSolver>optionValuePairs">List of option/value pairs</a>
     <a href="matlab: help helpPbeSolver>variables">Nomenclature  for input/output structures</a>
     <a href="matlab: help helpPbeSolver>inputStruct">Input structure: pbeDefinition</a>
       > <a href="matlab: help helpPbeSolver>seedingStruct">Structure for a seeding event.</a>
     <a href="matlab: help helpPbeSolver>outputStruct">Output structure array: result</a>

 See also: extractDataFromPbeResult, executePbeCase, loadPbeREsult,
   limitGrowthRate

This function calls:

addNucleation add nucleation handling to solver structure
checkInputSolver check fields in input data that must exist.
ensureProperDecData check if unified partitions are right and correct of required
genericDerivative compute derivative of states from structures
genericDerivativeStateVector compute derivative of states
genericEventFunction evaluate if events have to be triggered
newBinForNucleation remove nucleation handling to solver structure
obtainAbsTolVector create abs tol vector from given solver state structure
obtainCleanSolverState obtain an empty solver state with appropriate filling
packState pack solver state variables in a single vector
prepareSeeding ensure proper seeding structure
removeNucleation remove nucleation handling to solver structure
unpackState unpack state vector according to packing parameters

This function is called by:

executePbeCase execute a case from the database

0001 function [result pbeDefinition] = solvePbe(pbeDefinition, varargin)
0002 % solve a population balance (growth, nucleation, seeding, bulk states)
0003 %
0004 % Syntax: result = solvePopulationBalance(pbeDefinition)
0005 %
0006 % This routine essentially implements the method of characteristics for a
0007 %   set of input points. In this scope, the routine considers the
0008 %   limitation of the growth rate to the gemetrical boundaries (e.g.
0009 %   validity). It does also recognize dissolving particles.
0010 % Not explicitely contained in the scheme of this function is the
0011 %   integration of the particle population to, e.g. the total particle
0012 %   volume. This functionality is implicitely provided by the definition of
0013 %   bulk properties in the PBE definition structure (pbeDef). However, this
0014 %   library follows the approach of Monte Carlo integral estimates which
0015 %   is, in particular, reflected by the function extractDataFromPbeResult.
0016 % The support for multiple seeding and nucleation was cancelled because the
0017 %   overal scheme with Monte Carlo integral estimates does not yet cover
0018 %   these situations. It would, however, be straightforward to add this
0019 %   functionality. Your turn! ;)
0020 % For input/output see the link below.
0021 %
0022 % Detailed documentation for inpu/output structures and option/value pairs
0023 %   for the function call are available, seperately, <a href="matlab: helpPbeSolver">here</a>.
0024 %   This documentation contains information on:
0025 %     <a href="matlab: help helpPbeSolver>optionValuePairs">List of option/value pairs</a>
0026 %     <a href="matlab: help helpPbeSolver>variables">Nomenclature  for input/output structures</a>
0027 %     <a href="matlab: help helpPbeSolver>inputStruct">Input structure: pbeDefinition</a>
0028 %       > <a href="matlab: help helpPbeSolver>seedingStruct">Structure for a seeding event.</a>
0029 %     <a href="matlab: help helpPbeSolver>outputStruct">Output structure array: result</a>
0030 %
0031 % See also: extractDataFromPbeResult, executePbeCase, loadPbeREsult,
0032 %   limitGrowthRate
0033 
0034 % The elk-library: convex geometry applied to crystallization modeling.
0035 %   Copyright (C) 2012 Alexander Reinhold
0036 %
0037 % This program is free software: you can redistribute it and/or modify it
0038 %   under the terms of the GNU General Public License as published by the
0039 %   Free Software Foundation, either version 3 of the License, or (at your
0040 %   option) any later version.
0041 %
0042 % This program is distributed in the hope that it will be useful, but
0043 %   WITHOUT ANY WARRANTY; without even the implied warranty of
0044 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0045 %   General Public License for more details.
0046 %
0047 % You should have received a copy of the GNU General Public License along
0048 %   with this program.  If not, see <http://www.gnu.org/licenses/>.
0049 
0050 % start timer
0051 tic
0052 
0053 %% Input handling
0054 % variable input stuff
0055 optionStruct = mapOptionStruct(varargin, ...
0056     'odeSolver', @ode45, ...
0057     'absTolPosition', 1e-12, ...
0058     'absTolNumber', 1e10, ...
0059     'absTolBulk', 1e-12, ...
0060     'absTolDensity', 1e8, ...
0061     'relTol', 1e-5, ...
0062     'initialTimeStep', 1, ...
0063     'distanceLimit', 1e-8, ...
0064     'binLoad', 1e17, ...
0065     'gradientRelativeDelta', 1e-3, ...
0066     'gradientAbsoluteDelta', 1e-10, ...
0067     'updateUnifiedLimit', 1, ...
0068     'updateTiming', 1, ...
0069     'limitGrowthRate', 1, ...
0070     'useEmbeddingPartitions', 1, ...
0071     'verbose', 2, ...
0072     'density', 0);
0073 
0074 % allow pbeDef to overwrite some of the options above
0075 if isfield(pbeDefinition, 'optionStruct')
0076     pbeDefinition.optionStruct = mergeStruct(pbeDefinition.optionStruct, ...
0077         optionStruct);
0078 else
0079     pbeDefinition.optionStruct = optionStruct;
0080 end
0081 optionStruct = pbeDefinition.optionStruct;
0082 
0083 % prepare additional information
0084 pbeDefinition = checkInputSolver(pbeDefinition);
0085 % pbeDefinition = computeLimitCone(pbeDefinition);
0086 [pbeDefinition.seeding seedingTimeVector] = prepareSeeding(...
0087     pbeDefinition.seeding, pbeDefinition, optionStruct);
0088 
0089 %% Init: event handling
0090 nextSeedIndex = 1;  % index to next seeding population
0091 eventIndex = [];    % index of an event to handle
0092 % set seeding time vector: last event is simulation end
0093 seedingTimeVector = [seedingTimeVector pbeDefinition.timeSpan(end)];
0094 %% Init: solver state
0095 % use a central definition of al properly sorted fields
0096 solverState = obtainCleanSolverState;
0097 % time
0098 solverState.time = pbeDefinition.timeSpan(1);
0099 solverState.evalTime = toc;
0100 % state variables
0101 solverState.bulkStateVector = pbeDefinition.bulkInitial(:);
0102 solverState.hcMatrix = zeros(pbeDefinition.nDecCoordinate + ...
0103                              pbeDefinition.nExtraCoordinate, 0);
0104 solverState.numberVector = zeros(1, 0);
0105 if optionStruct.density
0106     solverState.probabilityVector = zeros(1, 0);
0107 end
0108 % help variables
0109 solverState.indexVector = zeros(1,0);
0110 solverState.upVector = zeros(1, 0);
0111 solverState.dissolvedPivotVector = [];
0112 % property information
0113 solverState.pivotPropertyMatrix = [];
0114 solverState.bulkPropertyVector = [];
0115 % option/structure information
0116 solverState.flagNucleation = 0;
0117 solverState.flagDensity = optionStruct.density;
0118 
0119 %% Init: nucleation state
0120 % if nucleation should be swiched on or off corresponds to the third event
0121 %   of the genericEventFunction
0122 [initialEventVector, ~, ~] = genericEventFunction(packState(solverState), ...
0123     solverState.time, pbeDefinition, solverState);
0124 if initialEventVector(1) > 0
0125     solverState = addNucleation(solverState, pbeDefinition);
0126 elseif initialEventVector(1) == 0
0127     error('elk:populationBalance:initialNucleation', ...
0128         ['The initial nucleation rate is exactly zero. ' ...
0129          'Please alter the initial nucleation or time slightly to resolve ' ...
0130          'this issue. If you don''t have nucleation, choose a negative ' ...
0131          'nucleation rate']);
0132 end
0133 
0134 %% Init: result
0135 % initialize field names
0136 result = solverState;
0137 result.hcMatrixDynamic = result.hcMatrix;
0138 result.derivativeStruct = struct([]);
0139 % ensure no data is written
0140 result(:) = [];
0141 % extend to target result vector
0142 result(length(pbeDefinition.timeSpan)).time = [];
0143 % fill time with nan for later filtering
0144 [result(:).time] = deal(nan);
0145 
0146 %% solver loop: BEGIN
0147 % 1) Add seeding pivots (if required)
0148 % 2) Setup ODE-state vector:
0149 %    * initial condition
0150 %    * handle for right-hand side
0151 % 3) Integrate to: (1) next seeding event,
0152 %                  (2) first nucleation bin full // start nucleation
0153 %                  (3) dissolution occures (unknown in which population)
0154 %                  (4) update unified partition data
0155 %                  (5) end of simulation
0156 % 4) Save solver states
0157 
0158 if optionStruct.verbose,disp('Start solving pbe ..');end
0159 while solverState.time < seedingTimeVector(end)
0160        
0161     %% solver loop:event handling
0162     eventIndex = unique(eventIndex);
0163     % correct event indexing for dissolved pivot Indexing
0164     %                original
0165     %                (in event function)
0166     % nuc on/off     1
0167     % nuc new bin    2
0168     % update up      3
0169     % confined error 4
0170     % dissolution    5..
0171     particleDissolveIndex = eventIndex(eventIndex > 4) - 4;
0172     while ~isempty(eventIndex)
0173         
0174         %% solver loop:switch nucleation on/off
0175         if eventIndex(1) == 1
0176             if solverState.flagNucleation
0177                 disp('  ..stopping nucleation');
0178                 solverState = removeNucleation(solverState);
0179                 % new bin event would be superflous now:
0180                 eventIndex(eventIndex == 2) = [];
0181                 % if simultaneaosly pivots grow out then there new
0182                 %   corresponding index is (boundary removed):
0183                 particleDissolveIndex = particleDissolveIndex - 1;
0184             else
0185                 disp('  ..starting of nucleation');
0186                 solverState = addNucleation(solverState, pbeDefinition, ...
0187                     optionStruct);
0188                 % if simultaneaosly pivots grow out then there new
0189                 %   corresponding event index is (bin + boundary added):
0190                 particleDissolveIndex = particleDissolveIndex + 2;
0191             end
0192             
0193         %% solver loop:nucleation new bin event
0194         elseif eventIndex(1) == 2
0195             disp(['  ..new nucleation bin, now: ' num2str(sum(...
0196                 solverState.indexVector==-1)+2)]);
0197             solverState = newBinForNucleation(solverState, pbeDefinition, ...
0198                 optionStruct);
0199             % if simultaneaosly pivots grow out then there new
0200             %   corresponding event index is (pivot added):
0201             particleDissolveIndex = particleDissolveIndex + 1;
0202             
0203         %% update unified partitions
0204         elseif eventIndex(1) == 3
0205             disp('  ..updating unified partition data') 
0206             [dump solverState.upVector] = ensureProperDecData(...
0207                 solverState.hcMatrix, solverState.upVector, pbeDefinition);
0208             
0209         %% confined vector created
0210         elseif eventIndex(1) == 4
0211             %! DEV: Yes, this part is not implemented because it was not
0212             %       necessary for case studies. See the
0213             %       thesis/documentation for details. The idea is to handle
0214             %       confined vectors by the closest unconfined partition.
0215             %       This is not as trivial as it might sound.
0216             error('elk:pbeSolver:confinedVectorGenerated', ['This ' ...
0217                 'simulation generated an unconfined vector. This can ' ...
0218                 'happen but you might try an implicit solver or different ' ...
0219                 'initial conditions to overcome this issue. An automatic ' ...
0220                 'handling of this case is not yet implemented. Sorry.']);
0221             %/
0222             
0223             
0224         %% solver loop:dissolution event
0225         % nucleatio pivot and boundary are treated otherwise, see the end
0226         %   of genericEventFunction for more information
0227         elseif eventIndex > 4 && ...
0228                particleDissolveIndex >= (1+2*solverState.flagNucleation) && ...
0229                particleDissolveIndex <= length(solverState.numberVector)
0230             % remove pivot position
0231             disp(['  ..removing pivot' num2str(particleDissolveIndex(1)) ...
0232                   ' .. remaining ' ...
0233                   num2str(length(solverState.numberVector)-1)]);
0234             % increase count for removed pivots
0235             thisPivotClass = solverState.indexVector(particleDissolveIndex(1));
0236             if thisPivotClass > 0
0237                 solverState.dissolvedPivotVector(thisPivotClass) = ...
0238                       solverState.dissolvedPivotVector(thisPivotClass) + 1;
0239             end
0240             solverState.indexVector(particleDissolveIndex(1)) = [];
0241             solverState.hcMatrix(:, particleDissolveIndex(1)) = [];
0242             solverState.upVector(particleDissolveIndex(1)) = [];
0243             % remove pivot numbers
0244             solverState.numberVector(particleDissolveIndex(1)) = [];
0245             if solverState.flagDensity
0246                 solverState.probabilityVector(particleDissolveIndex(1)) = [];
0247             end
0248             % we start with the first pivot so that ALL remaining pivots
0249             %   are indexed now as:
0250             particleDissolveIndex = particleDissolveIndex - 1;
0251             % remove particle index that was lined up:
0252             particleDissolveIndex(1) = [];
0253         %% no valid event
0254         else
0255             error('elk:solver:unknownEvent', 'This should not happen.');
0256         end
0257         %% remove event
0258         eventIndex(1) = [];
0259     end
0260 
0261     %% solver loop: seeding event
0262     if solverState.time >= seedingTimeVector(nextSeedIndex)
0263         % shortcuts
0264         thisHcMatrix = pbeDefinition.seeding(nextSeedIndex).positionMatrix;
0265         thisNumberVector = pbeDefinition.seeding(nextSeedIndex).numberVector;
0266         nDecCoordinate = pbeDefinition.nDecCoordinate;
0267         
0268         % filter useful points (having some particle number)
0269         newPivotFilter = thisNumberVector ~= 0;
0270         nNew = sum(newPivotFilter);
0271         nVoid = sum(~newPivotFilter);
0272         
0273         % assign hcMatrix
0274         solverState.hcMatrix(:, end + (1:nNew)) = thisHcMatrix(:, ...
0275             newPivotFilter);
0276         % assign particle number
0277         solverState.numberVector(1, end + (1:nNew)) = ...
0278             pbeDefinition.seeding(nextSeedIndex).numberVector(newPivotFilter);
0279         
0280         % assign upVecctor
0281         [hcMatrixDump solverState.upVector(1, end + (1:nNew)) ~] = ensureProperDecData(...
0282              thisHcMatrix(:, newPivotFilter), ones(1, nNew), pbeDefinition);
0283         
0284         %! Dev: below is obsolete after adding the above
0285 %         if pbeDefinition.dec.isComplete
0286 %             pivotPositionMeasureMatrix = mapToConfinementCone(...
0287 %                 thisHcMatrix(1:nDecCoordinate, newPivotFilter), ...
0288 %                 pbeDefinition.dec, 0);
0289 %             solverState.upVector(1, end + (1:nNew)) = computeUpVector(...
0290 %                 pbeDefinition, pivotPositionMeasureMatrix, ...
0291 %                 elkZerotol);
0292 %         else
0293 %             pivotPositionMeasureMatrix = ...
0294 %                 thisHcMatrix(1:nDecCoordinate, newPivotFilter);
0295 %             solverState.upVector(1, end + (1:nNew)) = 1;
0296 %         end
0297         %/
0298         
0299 
0300         solverState.indexVector(1, end + (1:nNew)) = ...
0301             ones(1, nNew)*nextSeedIndex;
0302         
0303         solverState.dissolvedPivotVector(end+1) = nVoid;
0304         
0305         if solverState.flagDensity
0306             solverState.probabilityVector(1, end + (1:nNew)) = ...
0307                 pbeDefinition.seeding(nextSeedIndex).probabilityVector(...
0308                   newPivotFilter);
0309         end
0310         
0311         nextSeedIndex = nextSeedIndex + 1;
0312     end
0313     
0314 %     disp('DBG: ...')
0315 %     solverState
0316 %     derivativeState = genericDerivative(solverState, 0, pbeDefinition)
0317 %     disp('..DBG/')
0318     
0319     %% solve a step
0320     thisStateVector = packState(solverState);
0321     % break if state vector is empty
0322     if isempty(thisStateVector)
0323         warning('elk:solver:emptyState', ...
0324             'Solver cancelled calculations because of empty state vector');
0325         return
0326     end
0327     
0328     thisDerivativeFunction = @(time, state)(genericDerivativeStateVector(...
0329         state, time, pbeDefinition, solverState));
0330 
0331     thisEventFunction = @(time, state)(genericEventFunction(state, time, ...
0332         pbeDefinition, solverState));
0333     
0334     thisTimeSpan = [solverState.time seedingTimeVector(nextSeedIndex)];
0335     % limit maximum step to updateTiming from total interval
0336     thisTimeSpan(2) = min(thisTimeSpan(2), thisTimeSpan(1) + ...
0337         optionStruct.updateTiming*...
0338         (seedingTimeVector(end) - seedingTimeVector(1)));
0339     thisTimeSpan = [thisTimeSpan(1), ...
0340                     pbeDefinition.timeSpan((pbeDefinition.timeSpan > thisTimeSpan(1)) & ...
0341                                            (pbeDefinition.timeSpan < thisTimeSpan(2))), ...
0342                     thisTimeSpan(2)];
0343     
0344     thisAbsTol = obtainAbsTolVector(solverState, ...
0345         optionStruct.absTolPosition, ...
0346         optionStruct.absTolNumber, ...
0347         optionStruct.absTolBulk, ...
0348         optionStruct.absTolDensity);
0349     
0350     
0351     [thisTimeVector, thisStateMatrix, eventTime, eventState, eventIndex] = ...
0352         feval(optionStruct.odeSolver,thisDerivativeFunction, thisTimeSpan, thisStateVector, ...
0353           odeset(...
0354             'absTol', thisAbsTol, ...
0355             'relTol', optionStruct.relTol, ...
0356             'events', thisEventFunction, ...
0357             'InitialStep', optionStruct.initialTimeStep ...
0358           ));
0359     if ~any([1 eventIndex] == 3)
0360         eventIndex = [eventIndex 3];
0361     end
0362 
0363     solverState = unpackState(thisStateMatrix(end, :)', thisTimeVector(end), ...
0364         solverState);
0365     if isempty(eventIndex)
0366         disp(['..stopped ode45 at time ' num2str(solverState.time)])
0367     else
0368         disp(['..stopped ode45 for event at time ' num2str(solverState.time) ...
0369             ' with event index ' num2str(eventIndex(:)')])
0370     end
0371     
0372     %% save unpacked states
0373     for iTime = 1:length(thisTimeVector)
0374         % hide event times where the solver also stops
0375         thisWriteIndex = find(thisTimeVector(iTime) == pbeDefinition.timeSpan, 1);
0376         if isempty(thisWriteIndex)
0377             % this is no time to be written
0378             continue
0379         end
0380         
0381         % initialize thisResult
0382         thisResult = unpackState(thisStateMatrix(iTime, :)', thisTimeVector(iTime), ...
0383             solverState);
0384         % save state values of hc-matrix
0385         thisResult.hcMatrixDynamic = thisResult.hcMatrix;
0386         % try to use previous unified index vector for speedup
0387         if thisWriteIndex > 1 && ~isnan(result(thisWriteIndex - 1).time)
0388             thisResult.upVector = result(thisWriteIndex - 1).upVector;
0389         end
0390         % compute derivative before it is too late (hcMatrix changed below)
0391         thisResult.derivativeStruct = genericDerivative(thisResult, ...
0392             thisTimeVector(iTime), pbeDefinition);
0393 
0394         % ensure proper data
0395         [thisResult.hcMatrix thisResult.upVector] = ...
0396             ensureProperDecData(thisResult.hcMatrix, thisResult.upVector, ...
0397             pbeDefinition);
0398         
0399         % pivot properties
0400         thisResult.pivotPropertyMatrix = pbeDefinition.pivotProperty(...
0401             thisResult.hcMatrix, thisResult.upVector);
0402         
0403         % bulk properties
0404         thisResult.bulkPropertyVector = pbeDefinition.bulkProperty(...
0405             thisResult.pivotPropertyMatrix, ...
0406             thisResult.numberVector, thisResult.bulkStateVector, ...
0407             thisResult.time);
0408         
0409         % assign
0410         thisResult.evalTime = toc;
0411         result(thisWriteIndex) = thisResult;
0412     end
0413     
0414 end
0415 disp('integration complete')
0416 
0417 %% cleanup results
0418 removeFilter = isnan([result(:).time]);
0419 if any(removeFilter)
0420     warning('elk:pbeSolver:missingOutput', ['Not all requested output '...
0421         'times were generated. Check [result(:).time] for ' ...
0422         num2str(sum(removeFilter)) ' missing entries.'])
0423 end
0424 result(removeFilter) = [];

solvePbe

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE