From 637112c702f83feb36936040ccbd4bd9ace2b63c Mon Sep 17 00:00:00 2001 From: VvanderHeiden Date: Tue, 27 Aug 2024 10:42:29 +0200 Subject: [PATCH 1/2] Reselove missing 'gamms_P' in workers bug (#53) --- bin/panther.m | 3 +-- examples/Mmax/example_Groningen_fault_stresses.m | 5 +++-- examples/example_depth_dependent_property.m | 4 ++-- src/calculator/FaultStress.m | 3 +++ src/calculator/FaultStressChange.m | 2 +- 5 files changed, 10 insertions(+), 7 deletions(-) diff --git a/bin/panther.m b/bin/panther.m index 6ad61f4..ecde24a 100644 --- a/bin/panther.m +++ b/bin/panther.m @@ -78,7 +78,7 @@ % of seconds. If it is already running (type gcp to check), parfor will initiate much faster parfor (i = 1 : n_members, matlab_workers) % disp([num2str(i),'/', num2str(n_members)]); - L{i} = y./sin(ensemble{i}.dip*pi/180); + L{i} = y./sin(ensemble{i}.dip*pi/180); % initial stress initial_stress{i} = InitialStress(y, ensemble{i}); @@ -126,7 +126,6 @@ reduced_load_table = analysis.load_table(indices_for_saving,:); - for i = 1 : n_members % save stresses, pressure, slip diff --git a/examples/Mmax/example_Groningen_fault_stresses.m b/examples/Mmax/example_Groningen_fault_stresses.m index 17dd16d..0c7849e 100644 --- a/examples/Mmax/example_Groningen_fault_stresses.m +++ b/examples/Mmax/example_Groningen_fault_stresses.m @@ -2,8 +2,9 @@ % TODO make into a function which uses a config file %% prepare input -proj = matlab.project.rootProject; -proj_folder = convertStringsToChars(proj.RootFolder); +% proj = matlab.project.rootProject; +% proj_folder = convertStringsToChars(proj.RootFolder); +proj_folder = '/Users/4146751/surfdrive2/PhD/DeepNL_Phys_Mmax/WP1/Code/PhysMmax_new/PANTHER'; % out_folder = [proj_folder,'/results/test_single_fault_segment/']; file_name.faults = [proj_folder,'/examples/Mmax/Mmax_files/Bourne_2017_FaultModel_Geometries_DuplicatesRemoved.xlsx']; diff --git a/examples/example_depth_dependent_property.m b/examples/example_depth_dependent_property.m index 7dd4637..64a92f6 100644 --- a/examples/example_depth_dependent_property.m +++ b/examples/example_depth_dependent_property.m @@ -77,8 +77,8 @@ % set depth varying dip sc2_variable_dip.input_parameters.dip.uniform_with_depth = 0; % make dip depth-variable y = sc2_variable_dip.y; -sc2_variable_dip.input_parameters.dip.value_with_depth = 90*ones(size(y)); -i_mid = floor(length(y)/2); +sc2_variable_dip.input_parameters.dip.value_with_depth = 80*ones(size(y)); +i_mid = floor(length(y)/2) sc2_variable_dip.input_parameters.dip.value_with_depth(i_mid:end) = 60; % sc2_variable_dip.input_parameters.dip.value_with_depth = 60 + y*dip_gradient_per_m; diff --git a/src/calculator/FaultStress.m b/src/calculator/FaultStress.m index 3084a73..23c062d 100644 --- a/src/calculator/FaultStress.m +++ b/src/calculator/FaultStress.m @@ -30,6 +30,9 @@ self.tau = initial.tau0 + change.dtau; end + function self = get_reactivation_nucleation_stress(self, reactivation_load_step, nucleation_load_step) + [self.sne_reac, self.tau_reac] = self.get_stress_at_load_step(reactivation_load_step); + end function self = get_nucleation_stress(self, nucleation_load_step) % INPUT diff --git a/src/calculator/FaultStressChange.m b/src/calculator/FaultStressChange.m index 27f6738..78cf923 100644 --- a/src/calculator/FaultStressChange.m +++ b/src/calculator/FaultStressChange.m @@ -123,7 +123,7 @@ end if length(gamma_PT) == 1 gamma = gamma_PT; - elseif length(gamma_P) == length(y) + elseif length(gamma_PT) == length(y) gamma = gamma_PT(j); end [dsn_j, dtau_j] = self.get_stress_change_component(GF{j}, dPT_FW, dPT_HW, gamma); From 14bc8497a03ef3717fb35b85bbdfe9c76e32ab72 Mon Sep 17 00:00:00 2001 From: VvanderHeiden Date: Tue, 27 Aug 2024 16:01:37 +0200 Subject: [PATCH 2/2] fixing heterogeneous elastic paramterization issues related to mu_II --- bin/MultiFaultCalculator.m | 58 +++++++++++++++++++++++++++----------- src/calculator/FaultSlip.m | 16 ++++++++--- 2 files changed, 53 insertions(+), 21 deletions(-) diff --git a/bin/MultiFaultCalculator.m b/bin/MultiFaultCalculator.m index fab58be..53da9ce 100644 --- a/bin/MultiFaultCalculator.m +++ b/bin/MultiFaultCalculator.m @@ -115,38 +115,62 @@ end function self = set_run_setting(self, setting_name, setting_value) - % specify run settings per pillar. (exc load table, y, input) + % Specify run settings per pillar. % INPUT - % setting name - % setting value cell array, or array of floats or single - % cell/single float + % setting_name - Name of the setting to be applied + % setting_value - Cell array, array of floats, single cell, single float, or string + [is_valid, value_type] = self.is_valid_setting_name(setting_name); + if is_valid + valid_value = 1; % Initialize as valid + for i = 1 : length(self.pillars) - valid_value = 1; - if iscell(setting_value) && (length(setting_value) == length(self.pillars)) - assign_value = setting_value{i}; - elseif iscell(setting_value) && (length(setting_value) == 1) - assign_value = setting_value{1}; - elseif isfloat(setting_value) && (length(setting_value) == 1) - assign_value = setting_value(1); - elseif isfloat(setting_value) && (length(setting_value) == self.n_pillars) - assign_value = setting_value(i); + if iscell(setting_value) + % Case 1: Cell array with the same length as pillars + if length(setting_value) == length(self.pillars) + assign_value = setting_value{i}; + % Case 2: Single cell element (could be string, char array, or numeric) + elseif length(setting_value) == 1 + assign_value = setting_value{1}; + else + valid_value = 0; % Invalid case: cell array with incorrect length + end + + elseif isnumeric(setting_value) + % Case 3: Single numeric value + if length(setting_value) == 1 + assign_value = setting_value(1); + % Case 4: Numeric array with the same length as pillars + elseif length(setting_value) == self.n_pillars + assign_value = setting_value(i); + else + valid_value = 0; % Invalid case: numeric array with incorrect length + end + + elseif ischar(setting_value) || isstring(setting_value) + % Case 5: String or character array + assign_value = setting_value; + else - valid_value = 0; + valid_value = 0; % Invalid case: unsupported type end - if valid_value & strcmp(value_type,class(assign_value)) + + % Assign value if valid and type matches the expected type + if valid_value && strcmp(value_type, class(assign_value)) self.pillars{i}.(setting_name) = assign_value; end end - if ~valid_value || ~strcmp(value_type,class(assign_value)) + + % Display a message if the value was not valid or the type didn't match + if ~valid_value || ~strcmp(value_type, class(assign_value)) disp('Specified setting type does not seem the right type or dimension, check'); end end - end + function [nuc_load_step] = get_minimum_nucleation_load_step(self) if self.run_done nuc_load_step = min(self.result_summary.nucleation_load_step); diff --git a/src/calculator/FaultSlip.m b/src/calculator/FaultSlip.m index b5a1601..1ca62fc 100644 --- a/src/calculator/FaultSlip.m +++ b/src/calculator/FaultSlip.m @@ -98,7 +98,7 @@ % length(time). slip_zone_indices = self.get_slip_zone_indices(slipping); slip_zone_length = nan(size(slip_zone_indices)); % along-fault length of slip zones - nucleation_length = nan(size(slip_zone_indices)); % nucleation length, per slip zone + nucleation_length = nan(size(slip_zone_indices)); % nucleation length, per slip zone % iterate over time steps for i = 1 : size(slip_zone_indices,2) % iterate over number of slip zones @@ -115,6 +115,7 @@ slip_zone_length(j,i) = L_start - L_end; sne_slip = (sne(slip_zone_indices{j,i}, i)); tau_slip = (tau(slip_zone_indices{j,i}, i)); + if length(f_s) == size(tau,1) % add ,i here if friction changes per timestep f_s_slip = f_s(slip_zone_indices{j,i}) ; @@ -131,9 +132,15 @@ else d_c_slip = d_c; end + if length(mu_II) == size(tau,1) + mu_II_slip = mean(mu_II(slip_zone_indices{j,i})); + else + mu_II_slip = mu_II; + end delta_tau = mean(sne_slip .* (f_s_slip - f_d_slip)); tau_0_d = mean(tau_slip - sne_slip .* f_d_slip); - nucleation_length(j,i) = self.calculate_nucleation_length( delta_tau, tau_0_d, d_c_slip, mu_II, nuc_crit, nuc_len_fixed); + + nucleation_length(j,i) = self.calculate_nucleation_length(delta_tau, tau_0_d, d_c_slip, mu_II_slip, nuc_crit, nuc_len_fixed); %nucleation_length(j,i) = 1.158 * mu_II * d_c./((f_s - f_d) * average_sne_in_slip_zone); end end @@ -210,8 +217,9 @@ % mu_II mode II shear stiffness nx = length(L); Kline = -nx/(2*pi*(L(1)-L(end))) ./ ( [0:nx-1]'.^2-0.25 ); - Ko = toeplitz(Kline); - K = mu_II*Ko; + Ko = toeplitz(Kline); + % K = mu_II.*Ko; + K = Ko*mu_II; K(and(K<0,K>(min(min(K)/10000)))) = 0; % set very small changes to 0 to avoid continued interactions of the two peaks to end