Merge pull request #84 from TNO/7-validation-of-input

7 validation of input

bin/panther.m

@@ -20,6 +20,9 @@
         analysis.generate_ensemble;
     end
 
+    % validate input
+    validate_input(analysis);
+
     % unpack some parameters
     y = analysis.y;
     load_table = analysis.load_table;
@@ -32,8 +35,8 @@
     ensemble = analysis.ensemble;
     nucleation_criterion = analysis.nucleation_criterion;
     nucleation_length = analysis.nucleation_length_fixed;
-    dy = y(1) - y(2);
- 
+
+
     % define output steps
     if ~ismember(analysis.save_stress,'none')
            if ismember(analysis.save_stress,'all')

resources/project/R6m7AxFfSX7StHtVJ9WG44h96GM/xMbsyEyRLQ9USX_pduFdZqREeckd.xml

@@ -0,0 +1,6 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<Info>
+    <Category UUID="FileClassCategory">
+        <Label UUID="design"/>
+    </Category>
+</Info>

resources/project/R6m7AxFfSX7StHtVJ9WG44h96GM/xMbsyEyRLQ9USX_pduFdZqREeckp.xml

@@ -0,0 +1,2 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<Info location="validate_input.m" type="File"/>

src/calculator/Greens_functions_stress/Gxx_combined.m

@@ -0,0 +1,80 @@
+function Gxx = Gxx_combined(a,b,c,d,theta,x,y)
+    % Gxx_combined combines Green's functions for different reservoir compartments.
+    %
+    % This function combines Green's functions for different reservoir compartments.
+    % The reservoir is modeled through a combination of a rectangles and triangles.
+    %
+    % Parameters:
+    %   a - Vertical distance to inner reservoir corners, from center (y = 0)
+    %   b - Vertical distance to outer reservoir corners, from center (y = 0)
+    %   c - Length of the left side of the reservoir
+    %   d - Length of the right side of the reservoir 
+    %   theta - Angle of the fault [radians]
+    %   x - x-coordinate of evaluation point(s)
+    %   y - y-coordinate of evaluation point(s)
+    %
+    % Returns:
+    %   Gxx - Combined Green's function value
+    %
+    % Notes:
+    %   - If c = 0, the reservoir is on the right side only.
+    %   - If d = 0, the reservoir is on the left side only.
+    %   - A small number is added to x and y to avoid singularities and division by zero.
+    %
+    % Example:
+    %   Gxx = Gxx_combined(10, 20, 5, 5, pi/4, 2, 3);
+
+    % Combines Green's functions for the different reservoir compartments. 
+    % The reservoir is modeled through a combination of a rectangle and triangle 
+    % The function can be used for reservoirs on two sides of the fault or 
+    % one side
+    % c = 0 --> right side only
+    % d = 0 --> left side only
+
+    num_offset = 1e-9; % small number to add to x and y to avoid singularities and division by 0) [-]  
+    if x == 0 || x == d
+        x = x + num_offset*abs(d);
+        x = x + num_offset;
+    end
+    if x == -c 
+        x = x - num_offset*abs(c);
+        x = x - num_offset;
+    end
+    if y == a || y == b
+        y = y + num_offset*abs(b-a);
+    end
+    if y == -a || y == -b
+        y = y - num_offset*abs(b-a);
+    end
+
+    % Compute integrals of individual components
+    f = b/tan(theta);       % offset in x 
+    e = a/tan(theta);
+
+    if and(~(c==0), ~d==0)
+        Gxx_left_rectangle = Gxx_rectangle(-c,-f,-b,a,x,y);
+        Gxx_left_triangle  = Gxx_triangle(-e,f,-a,b,-x,-y);
+        Gxx_right_triangle = Gxx_triangle(-e,f,-a,b,x,y);
+        Gxx_right_rectangle = Gxx_rectangle(f, d,-a,b,x,y);
+    elseif c==0 && ~d==0
+        % right side only
+        Gxx_left_rectangle = 0;
+        Gxx_left_triangle  = 0;
+        Gxx_right_triangle = Gxx_triangle(-e,f,-a,b,x,y);
+        Gxx_right_rectangle = Gxx_rectangle(f, d,-a,b,x,y);
+    elseif ~c==0 && d==0
+         % left side only
+        Gxx_left_rectangle = Gxx_rectangle(-c,-f,-b,a,x,y);
+        Gxx_left_triangle  = Gxx_triangle(-e,f,-a,b,-x,-y);
+        Gxx_right_triangle = 0;
+        Gxx_right_rectangle = 0;
+    end
+
+    Gxx = Gxx_left_rectangle + Gxx_left_triangle + Gxx_right_rectangle + Gxx_right_triangle;
+
+end
+
+
+
+
+

src/calculator/ModelGeometry.m

@@ -30,10 +30,10 @@
     properties
         depth_mid (1,1) double {mustBeNegative} = -3000
         thick (1,1) double {mustBePositive} = 150
-        dip (:,1) double = 66
-        dip_azi (:,1) double = 0
+        dip (:,1) double {mustBeInRange(dip, 0, 90)} = 60
+        dip_azi (:,1) double = 90
         throw (1,1) double = 75
-        width_HW (1,1) = Inf
+        width_HW (1,1) {mustBeNonnegative} = Inf
         width_FW (1,1) = Inf
     end
 

src/calculator/PantherInput.m

@@ -11,10 +11,10 @@
         stochastic logical = 0;                     % activate stochastic analysis
         n_stochastic {mustBeInteger} = 1;           % number of stochastic runs
         save_stress cell = {'all'};                 % indicate which stress to save. 'all', 'none', 'first','last',[step_numbers]
-        load_case = 'P';                            % load case 'P': pressure changes, 'T': temperature changes
+        load_case {mustBeMember(load_case, {'P','T','PT'})} = 'P';               % load case 'P': pressure changes, 'T': temperature changes
         load_table table                            % table containing time steps, P and T steps (len(y), len(timesteps) for both FW and HW
         ensemble_generated = 0;                     % toggle specifying whether model ensemble has been generated
-        ensemble                                    % ensemble of n_stochastic members 
+        ensemble
         parallel logical = 1                        % parallel computing for large number of simulations
         aseismic_slip logical = 1                   % compute aseismic slip during nucleation phase
         nucleation_criterion {mustBeMember(nucleation_criterion,{'fixed','UR2D','Day3D','Ruan3D'})} = 'UR2D';   
@@ -25,6 +25,9 @@
         dx double  = 0;                             % [m] distance from from (for now only on fault allowed)
     end
 
+    properties (Dependent) 
+    end
+
     methods
 
         function self = PantherInput()
@@ -52,33 +55,6 @@ function reset_ensemble(self)
             self.ensemble_generated = 0;
         end
 
-        function self = generate_ensemble_from_table(self, input_table)
-            % generates a model ensemble from a table
-            % turn off stochastic (for now). Input values not in the input
-            % table will have a fixed value. No-depth dependency. 
-            self.stochastic = 0; 
-            % check for matching table field names
-            panther_input_names = properties(self.input_parameters);
-            table_names = fields(input_table);
-            matching_columns = ismember(table_names, panther_input_names);
-            n_ensemble = height(input_table);
-            column_indices = [];
-            if any(matching_columns)
-                column_indices = find(matching_columns);
-                self.ensemble_generated = 1;
-                for i = 1 : n_ensemble
-
-                    self.ensemble{i,1} = PantherMember(self.input_parameters, 0);
-                    for j = 1 : length(column_indices)
-                        var = table_names{column_indices(j)};
-                        self.ensemble{i,1}.(var) = input_table.(var)(i);
-                    end
-                end
-            else
-                disp('Table column names do not match Panther input parameter names');
-            end     
-        end
-
         function ensemble_table = ensemble_to_table(self)
             % create table of input parameter values
             props = properties(self.input_parameters);

src/calculator/PantherMember.m

@@ -4,28 +4,28 @@
     % singel number or an array of length(y)
 
     properties
-        young               % [Pa] Young's modulus
-        poisson             % [-] Poisson's ratio    
-        biot                % [-] Biot coefficient
-        therm_exp           % [1/K] thermal expansion coefficient
-        sH_dir              % [deg] direction SHmax from north
-        sHsh                % [-] ratio SH/Sh
-        shsv                % [-] ratio Sh/Sv
-        sv_grad             % [MPa/km] Vertical stress gradient
-        sv_offset           % [MPa] Offset vertical stress gradient at y=0
-        p_grad              % [MPa/km] pressure gradient
-        p_offset            % [MPa] offset of pressure gradient at y=0    
-        p_over              % [MPa] overpressure in the reservoir
-        p_grad_res          % [MPa/km] pressure gradient in reservoir  
-        hyd_diffusivity     % [m2/s] hydraulic diffusivity
-        T_grad              % [K/km] temperature gradient
-        T_offset            % [k] offset temperature gradient at y=0
-        dT_dy_multiplier    % [deg/m] multiplier dT in reservoir wr.t. reservoir mid. -ve is increasing with depth
-        therm_diffusivity 
-        f_s
-        f_d
-        d_c
-        cohesion
+        young double {mustBePositive} = 15e3                    % [MPa] Young's modulus
+        poisson double {mustBeInRange(poisson, 0, 0.5)} = 0.2   % [-] Poisson's ratio    
+        biot double {mustBeInRange(biot, 0, 1)} = 1             % [-] Biot coefficient
+        therm_exp double {mustBePositive} = 1e-5                % [1/K] thermal expansion coefficient
+        sH_dir double = 0                                       % [deg] direction SHmax from north
+        sHsh double {mustBePositive} = 1                        % [-] ratio SH/Sh
+        shsv double {mustBePositive} = 0.75                     % [-] ratio Sh/Sv
+        sv_grad double {mustBePositive} = 22                    % [MPa/km] Vertical stress gradient
+        sv_offset double = 0                                    % [MPa] Offset vertical stress gradient at y=0
+        p_grad double {mustBePositive} = 10.5                   % [MPa/km] pressure gradient
+        p_offset double = 0                                     % [MPa] offset of pressure gradient at y=0    
+        p_over double = 0                                       % [MPa] overpressure in the reservoir
+        p_grad_res double {mustBePositive} = 10.5               % [MPa/km] pressure gradient in reservoir  
+        hyd_diffusivity double {mustBePositive} = 1e-6          % [m2/s] hydraulic diffusivity
+        T_grad double {mustBePositive} = 31                     % [K/km] temperature gradient
+        T_offset double {mustBePositive} = 10                   % [k] offset temperature gradient at y=0
+        dT_dy_multiplier = 1                                    % [deg/m] multiplier dT in reservoir wr.t. reservoir mid. -ve is increasing with depth
+        therm_diffusivity double {mustBePositive} = 1e-6        % thermal diffusivity 
+        f_s double {mustBeInRange(f_s, 0, 1)} = 0.6             % [-] static friction coefficient
+        f_d double {mustBeInRange(f_d, 0, 1)} = 0.45            % [-] dynamic friction coefficient
+        d_c double {mustBePositive}                             % [m] critical slip distance for linear slip weakening function
+        cohesion double {mustBeNonnegative} = 0                 % [MPa] cohesion
     end    
 
     methods

src/calculator/PantherParameterList.m

@@ -11,7 +11,7 @@
         throw = PantherParam(50, 'Thickness','m', 1, nan, 0, 'uniform', 25, 75);                         % [m] vertical fault offset
         width_FW = PantherParam(inf, 'Width of footwall block','m', 1, nan, 0, 'uniform', 1000,5000 );   % [m] width footwall compartment
         width_HW = PantherParam(inf, 'Width of hanging wall block','m', 1, nan,0, 'uniform', 1000,5000);% [m]width footwall compartment
-        young = PantherParam(15e3, 'Young''s modulus','Pa', 1, nan, 0, 'uniform', 10e3, 20e3);           % [MPa] Young's modulus
+        young = PantherParam(15e3, 'Young''s modulus','MPa', 1, nan, 0, 'uniform', 10e3, 20e3);           % [MPa] Young's modulus
         poisson = PantherParam(0.15, 'Poisson''s ratio','-',1, nan, 0, 'uniform', 0.1, 0.2);             % [-] Poisson's ratio    
         biot = PantherParam(1, 'Biot coefficient','-', 1, nan, 0, 'uniform',0.7, 1.0);                   % [-] Biot coefficient
         therm_exp = PantherParam(1e-5, 'Thermal expansion coefficient','1/K', 1, nan, 0, 'uniform', 0.9e-5, 1.2e-5);  %[1/K] thermal expansion coefficient
@@ -33,6 +33,23 @@
         f_d =  PantherParam(0.45, 'Dynamic friction coefficient','-', 1, nan, 0, 'uniform',0.35,0.49);    % [-] dynamic friction coefficient
         d_c =  PantherParam(0.005, 'Critical slip distance','m', 1, nan, 0, 'uniform',  0.002, 0.01);     % [-] critical slip distance
         cohesion = PantherParam(0, 'Cohesion','MPa', 1, nan, 0, 'uniform', 0, 5);                         % [MPa] cohesion
-    end     
+    end
+
+    methods
+
+        function self = PantherParameterList()
+        end
+
+        function [depth_dependent_properties] = get_depth_dependent_properties(self)
+            props = properties(self);
+            depth_dependent_properties = {};
+            for i = 1 : length(props)
+                if self.(props{i}).uniform_with_depth == false
+                    depth_dependent_properties{end + 1} = props{i};
+                end
+            end
+        end
+
+    end
 
 end

src/validate_input.m

@@ -0,0 +1,94 @@
+function validate_input(analysis)
+    % function to validate input which has not yet been validated in the
+    % class properties, and perform sanity checks on the chosen input
+    % parameters
+
+    % check that the correct input is entered for the save_stress option
+    if ischar(analysis.save_stress{1})
+        if ~ismember(analysis.save_stress{1},{'all','none','first_last','first','last'});
+             error(['For save stress selection enter step numbers or ',...
+                'all, none, first_last, first, last' ]);
+        end
+    end
+
+    % check sensibility of the chosen Young's modulus
+    if log10(analysis.input_parameters.young.value) < 3 || log10(analysis.input_parameters.young.value) > 5
+        warning('NB default input unit Youngs modulus is MPa, please double check your input');
+    end
+
+    % check sensibility of the chosen vertical stress gradient
+    if analysis.input_parameters.sv_grad.value < 18 || analysis.input_parameters.sv_grad.value > 28
+        warning(['Vertical stress gradient in in MPa/km, and lies typically ',...
+            'in the range of 19 - 24 MPa/km. Given input was ', ...
+            num2str(analysis.input_parameters.sv_grad.value),' please check']);
+    end
+
+    % check sensibility of the chosen pressure gradient
+    if analysis.input_parameters.p_grad.value < 9 || analysis.input_parameters.p_grad.value > 15
+        warning(['Pore stress gradient in in MPa/km, and lies typically ',...
+            'in the range of 10 - 12 MPa/km. Given input was ', ...
+            num2str(analysis.input_parameters.p_grad.value),' please check']);
+    end
+
+    % check sensibility of the chosen temperature gradient
+    if analysis.input_parameters.T_grad.value < 26 || analysis.input_parameters.T_grad.value > 40
+        warning(['Pore stress gradient in in deg/km, and lies typically ',...
+            'in the range of 26 - 35 deg/km far away from plate boundaries.',...
+            'Given input was ', ...
+            num2str(analysis.input_parameters.T_grad.value),' please check']);
+    end
+
+    % check sensibility of the chosen depth
+    if analysis.input_parameters.depth_mid.value < -5000 
+        warning(['Depth of subsurface activities is typically < 5000 m deep ',...
+            'Given input was ', ...
+            num2str(analysis.input_parameters.depth_mid.value),' please check']);
+    end
+
+    % check for error in throw
+    if analysis.input_parameters.throw.value > abs(analysis.input_parameters.depth_mid.value)
+        error(['Throw of ', num2str(analysis.input_parameters.throw.value),...
+            ' exceeds reservoir depth']);
+    end
+
+    % check for error in throw
+    if analysis.input_parameters.thick.value > abs(analysis.input_parameters.depth_mid.value)
+        error(['Thickness of ', num2str(analysis.input_parameters.thick.value),...
+            ' exceeds reservoir depth']);
+    end
+    % add a time = 0 row to table if that is missing
+    if ~(analysis.load_table.time_steps(1) ==  0)
+        old_table = analysis.load_table;
+        first_row = analysis.load_table(1,:);
+        first_row{1,:} = 0;
+        analysis.load_table = [first_row; old_table];
+        disp('Added time step 0 row to load table');
+    end
+
+    depth_dependent_properties = analysis.input_parameters.get_depth_dependent_properties();
+    allowed_depth_dependent_properties = {'dip','dip_azi','young','poisson',...
+        'biot','therm_exp','sH_dir','shsv','sHsh','f_s','f_d','d_c','cohesion'};
+        if ~isempty(depth_dependent_properties)
+            faulty_input = ~ismember(depth_dependent_properties, allowed_depth_dependent_properties);
+            if any(faulty_input)
+                i_first_faulty = find(faulty_input, 1, 'first');
+                error(['Depth dependent property set where not allowed ',...
+                    'e.g. ', depth_dependent_properties{i_first_faulty}, ...
+                    '. Depth dependent properties are dip,dip_azi, young, poisson,',...
+                    'biot,therm_exp,sH_dir,shsv,sHsh,f_s, f_d,d_c,cohesion']);
+            end
+            % set stochastic mode for depth dependent properties off
+            for i = 1 : length(depth_dependent_properties)
+                input_value_length = length(analysis.input_parameters.(depth_dependent_properties{i}).value_with_depth);
+                if (length(analysis.y) ~= length(input_value_length)) & (~input_value_length == 1 )
+                    error(['Length of depth dependent property ', depth_dependent_properties{i},...
+                        ' does not seem to match length of y, please check input']);
+                elseif input_value_length == 1 
+                    warning(['Depth dependency set for ', depth_dependent_properties{i},...
+                        ' but length of value with depth is 1']);
+                end
+                analysis.input_parameters.(depth_dependent_properties{i}).stochastic = 0;
+            end
+        end
+
+end