evaluate if events have to be triggered Event list: 1) dissolution = distance to inner empty polytope cone 2) nucleation: new bin = binLoadMax - binLoad 3) nucleation: on/off = nucleation rate 4) confined vector = test for any(partition index = unconfined) THIS IS NO USER FUNCTION
0001 function [eventValue, eventTerminate, eventDirection] = ... 0002 genericEventFunction(state, time, pbeDefinition, solverStateOriginal) 0003 % evaluate if events have to be triggered 0004 % 0005 % Event list: 0006 % 1) dissolution = distance to inner empty polytope cone 0007 % 2) nucleation: new bin = binLoadMax - binLoad 0008 % 3) nucleation: on/off = nucleation rate 0009 % 4) confined vector = test for any(partition index = unconfined) 0010 % 0011 % THIS IS NO USER FUNCTION 0012 0013 % The elk-library: convex geometry applied to crystallization modeling. 0014 % Copyright (C) 2012 Alexander Reinhold 0015 % 0016 % This program is free software: you can redistribute it and/or modify it 0017 % under the terms of the GNU General Public License as published by the 0018 % Free Software Foundation, either version 3 of the License, or (at your 0019 % option) any later version. 0020 % 0021 % This program is distributed in the hope that it will be useful, but 0022 % WITHOUT ANY WARRANTY; without even the implied warranty of 0023 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0024 % General Public License for more details. 0025 % 0026 % You should have received a copy of the GNU General Public License along 0027 % with this program. If not, see <http://www.gnu.org/licenses/>. 0028 0029 %! Debug: 0030 %% Debug: visualize solver time stepping 0031 % % console output 0032 % disp(['Event.. ..with time: ' num2str(time, 16)]) 0033 % 0034 % % graphical output 0035 % persistent iCall lastTime 0036 % if isempty(iCall) 0037 % iCall = 0;lastTime = 0; 0038 % else 0039 % iCall = iCall + 1; 0040 % end 0041 % figure(1),hold on,plot(time, time-lastTime, 'b.', 'markerSize', 5) 0042 % lastTime = time; 0043 %/ 0044 0045 %% Info 0046 %! Dev: 0047 % Unfortunatelly, it is required to know the derivative. And to identify 0048 % disappearing facets, the unpacking is mandatory. 0049 % However, this function is planned for improvement. 0050 %/ 0051 nTotalCoordinates = pbeDefinition.nDecCoordinate + ... 0052 pbeDefinition.nExtraCoordinate; 0053 nDecCoordinates = pbeDefinition.nDecCoordinate; 0054 0055 %% event configuration 0056 eventTerminate = [ 1; 1; 1; 1; solverStateOriginal.numberVector'*0 + 1]; 0057 eventDirection = [ 0;-1; 0; 0; solverStateOriginal.numberVector'*0 + 1]; 0058 eventValue = nan(size(eventTerminate)); 0059 0060 %% generate derivative WITH nucleation activated (always) 0061 % to track the nucleation on/off event the nucleation must be switched on 0062 % at all times 0063 eventSolverState = unpackState(state, time, solverStateOriginal); 0064 if ~solverStateOriginal.flagNucleation 0065 0066 % in this case where nucleation is switched off, the pivot and boundary 0067 % position do NOT matter. The pivots and bins do exist only so that 0068 % the nucleation rate is calculated. Hence, they can simply be filled 0069 % by zeros. 0070 eventSolverState = addNucleation(eventSolverState, pbeDefinition); 0071 0072 % get derivative 0073 derivativeState = genericDerivative(eventSolverState, time, pbeDefinition); 0074 else 0075 % get derivative from original state/stuff 0076 derivativeState = genericDerivative(eventSolverState, time, pbeDefinition); 0077 end 0078 0079 %% start / stop of nucleation 0080 eventValue(1) = derivativeState.numberVector(1); 0081 0082 %% new nucleation bin 0083 if solverStateOriginal.flagNucleation 0084 eventValue(2) = pbeDefinition.optionStruct.binLoad ... 0085 - state(nTotalCoordinates + 1); 0086 else 0087 eventValue(2) = -1; 0088 end 0089 0090 %% update unified partition information 0091 eventValue(3) = derivativeState.fractionCorrected - ... 0092 pbeDefinition.optionStruct.updateUnifiedLimit; 0093 0094 %% confinement tracking 0095 if pbeDefinition.optionStruct.limitGrowthRate 0096 if any(pbeDefinition.dec.unifiedIsConfinedVector(... 0097 derivativeState.partitionIndexVector)) || ... 0098 any(derivativeState.partitionIndexVector == 0) 0099 eventValue(4) = 1; 0100 else 0101 eventValue(4) = -1; 0102 end 0103 else 0104 eventValue(4) = 1; 0105 end 0106 0107 0108 %% dissolution 0109 % If the vector grows out, the boundary cannot be mapped properly so that 0110 % the correction value from mapToValidCone is necessarily +Inf or -Inf. 0111 % The event simply jumps from -0.5 to 0.5 if the corresponding facet grows 0112 % out. 0113 if pbeDefinition.dec.isComplete 0114 emptyCone = obtainFacetValidityCone(pbeDefinition.dec, pbeDefinition.dec.nH+1); 0115 distanceFromEmptyConeVector = max(emptyCone.A * ... 0116 eventSolverState.hcMatrix(1:nDecCoordinates, :) + ... 0117 10*pbeDefinition.optionStruct.distanceLimit, [], 1); 0118 if ~isempty(pbeDefinition.extraZeroBound) 0119 % the distance from the extra condition must also be considered.. 0120 % this is an <the above> OR <the below> must be fulfilled for a 0121 % particle to remain existent 0122 distanceFromEmptyConeVector = min(distanceFromEmptyConeVector, ... 0123 max(pbeDefinition.extraZeroBound * ... 0124 eventSolverState.hcMatrix((nDecCoordinates+1):end), [], 1)); 0125 end 0126 else 0127 % there is no dissolution for incomplete decomposition data 0128 distanceFromEmptyConeVector = 1+0*eventSolverState.hcMatrix(1, :); 0129 end 0130 0131 if solverStateOriginal.flagNucleation 0132 eventValue(5:end) = distanceFromEmptyConeVector; 0133 elseif length(eventTerminate) > 3 0134 eventValue(5:end) = distanceFromEmptyConeVector(3:end); 0135 elseif length(eventTerminate) == 3 0136 % in this case, nucleation was switched off and no pivots exist so far 0137 else 0138 error('elk:populationBalance:internalError', 'This should not happen.'); 0139 end 0140 %! DEV: improve speed for dissolution events 0141 % eventValue(4:end) = sign(eventValue(4:end)).*abs(eventValue(4:end)).^0.1;