commit 5ecc6387c0125850ae904911cde191f9155a0e01 Author: Joanna Kocot Date: Fri Mar 31 13:00:31 2017 +0200 Stress Maps application added diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5881870 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +/* +/*/ +!.gitignore +!/src/ diff --git a/src/StressMaps/stress_maps_service.m b/src/StressMaps/stress_maps_service.m new file mode 100644 index 0000000..5290ba0 --- /dev/null +++ b/src/StressMaps/stress_maps_service.m @@ -0,0 +1,593 @@ +function stress_maps_service(root_folder, resolution, weightings, keyword, nparam, mc, stresstype, ldvalue, plot_fracs) + + % STRESS_MAPS_SERVICE Plot stress patterns associated with + % the hydraulic fracturing of a shale gas well. The data is based on + % the physical and geological conditions at the Preese Hall Episode site. + + % STRESS_MAPS_SERVICE(ROOT_FOLDER, RESOLUTION, WEIGHTINGS, KEYWORD, NPARAM, MC, STRESSTYPE, LDVALUE) + % plots a stress map at a high ([20]m) or low ([50]m) RESOLUTION for the STRESSTYPE + % ([N]ormal, [S]hear or [C]oulomb) using a given source WEIGHTING ratio. The parameter + % to be plotted is defined in the variable KEYWORD (ISIP, FlowRate, PumpTime) and the value + % of this in NPARAM. MC is the monte carlo iteration (1-50) and ROOT_FOLDER the location + % of all the mat files. The lateral distance is calculated using the + % LDVALUE threshold + + % Input options for nparam: + % dP = [0.426 0.676 0.926 1.176 1.676 1.926 2.176 2.426 2.926 3.426 3.676 4.426 4.926 5.176 6.676 8.176 9.676 11.176 12.676 14.176 15.676 17.176 18.676 20.176 21.676 23.176 24.676 26.176 27.676 29.176] + % FlowRate = [0.0983 0.1099 0.1215 0.1331 0.1446 0.1562 0.1678 0.1793 0.1909] + % PumpTime = [30 60 90 120 150 180 210 240 270 300 330 360] + + % version 1.1 20/12/2016 + % Copyright 2016 Rachel Westwood, Keele University + + % This work is licensed under a Creative Commons Attribution-ShareAlike + % 4.0 International License. For full details of the license see + % https://creativecommons.org/licenses/by-sa/4.0/legalcode + + + % Version History + % 1.1 - corrected colorbar issue with 20m resolution and filename issue + % with DP + + + % check the weightings to ensure they are correct + if length(weightings) ~= 3 + error('source weighting combination must contain 3 elements'); + end + + if sum(weightings) ~= 1 + error('source weighting combination must add up to 1'); + end + + if ldvalue < 0.001 || ldvalue > 1 + error('Lateral respect distance stress threshold should be between 0.001 and 1 MPa'); + end + + % check that nparam is in the list + dP = round(([37.25 37.5 37.75 38 38.5 38.75 39 39.25 39.75 40.25 40.5 41.25 41.75 42:1.5:66] - 36.824)*1000)/1000; % dP in MPa + flow = round((8500:1000:16500) * 1.157e-5 * 10000)/10000; % flow in m3/s + pumpt = 30:30:360; + + % if the keyword is differential pressure set it to ISIP. For dP and + % FlowRate adjust the nparam value. Check that the parameters are + % valid + nparam_orig = nparam; + keyword_orig = keyword; + switch upper(keyword) + case 'DP' + keyword = 'ISIP'; + if isempty(find(dP == nparam, 1)) + error('Parameter value invalid'); + end + nparam = (nparam + 36.824) * 1e6; + case 'FLOWRATE' + if isempty(find(flow == nparam, 1)) + error('Parameter value invalid'); + end + nparam = round((nparam * 86400)/100)*100; + case 'PUMPTIME' + if isempty(find(pumpt == nparam, 1)) + error('Parameter value invalid'); + end + end + + % generate empty matrices for the lateral distance calculation + all_x = []; + all_y = []; + all_z = []; + all_stress = []; + + % have to get the data for high and low res for the distance + % calculation, will only plot the required one. + + for interval = [20 50] + + %set the index + if interval == 20 + ind = 1; + else + ind = 2; + end + + % generate the filenames of each of the mechanisms for 100% i, 100% s + % and 100% d + f_i = fullfile(root_folder, ... + [keyword '_' num2str(nparam) '_i100-s0-d0_' ... + num2str(interval) 'm.mat']); + f_s = fullfile(root_folder, ... + [keyword '_' num2str(nparam) '_i0-s100-d0_' ... + num2str(interval) 'm.mat']); + f_d = fullfile(root_folder, ... + [keyword '_' num2str(nparam) '_i0-s0-d100_' ... + num2str(interval) 'm.mat']); + + % check the mat files are there + if ~exist(f_i, 'file') || ~exist(f_d, 'file') || ~exist(f_s, 'file') + error(['Not all mat files are present: ' num2str(nparam_orig)]); + end + + % Get the data from the mat files + m_i = load(f_i); + m_s = load(f_s); + m_d = load(f_d); + + % save the x, y, z + x_i = m_i.x; + y_i = m_i.y; + z_i = m_i.z; + + % x = x_i(:); + % y = y_i(:); + + + zz = zeros(length(x_i), 50); + + % get the correct stress matrix + switch (upper(stresstype)) + case 'C' + c_i = weightings(1) .* m_i.coulomb; + c_s = weightings(2) .* m_s.coulomb; + c_d = weightings(3) .* m_d.coulomb; + st = 'Coulomb'; + case 'N' + c_i = weightings(1) .* m_i.normal; + c_s = weightings(2) .* m_s.normal; + c_d = weightings(3) .* m_d.normal; + st = 'Normal'; + case 'S' + c_i = weightings(1) .* m_i.shear; + c_s = weightings(2) .* m_s.shear; + c_d = weightings(3) .* m_d.shear; + st = 'Shear'; + otherwise + error('Stress type not recognised'); + end + + % Add the stresses together + c = c_i + c_s + c_d; + + %preallocate the maximum stress matrix + lux = length(unique(x_i)); + luy = length(unique(y_i)); + stress_m = zeros(lux, luy, 50); + s_m = zeros(lux*luy, 50); + %xyz = zeros(50,1); + + % get the max stress for the required monte carlo iteration + for jj = 1:50 + [stress_m(:,:,jj), xyz(jj)] = max_stress(x_i, y_i, z_i, c(:,:,jj)); %#ok + s_m(:,jj) = reshape(stress_m(:,:,jj), [], 1); + zz(:,jj) = xyz(jj).z(:); + end + + all_x = [all_x; x_i(:)]; + all_y = [all_y; y_i(:)]; + all_z = [all_z; zz]; + all_stress = [all_stress; s_m]; + + % plot the stress map + if resolution == interval% get the fractures + fx = m_i.fracturesX; + fy = m_i.fracturesY; + fracX = fx{mc}; + fracY = fy{mc}; + fig = plot_stress(fracX, fracY, stress_m(:,:,mc), xyz(mc), plot_fracs, st, resolution); + end + + end + + % calculate the lateral distance + [mean_ld, median_d, max_ld, min_ld] = ... + lateralDistance(all_x, all_y, all_stress, ldvalue); + + % add the text with poisson ratio, youngs modulus, regional stress, + % flow, area and lateral distance + filename = [keyword_orig '_' st '_' num2str(resolution) 'm_' ... + num2str(nparam_orig) '_' ... + num2str(weightings(1)*100) 'i' num2str(weightings(2)*100) 's' num2str(weightings(3)*100) 'r__' ... + num2str(mc)]; + saveText(filename, m_i, keyword, [mean_ld, median_d, max_ld, min_ld], ldvalue, mc); + saveas(fig, [filename '.jpg']); +end + +function saveText(filename, m_i, keyword, latdist, latdistthreshold, mc) + % Function which generates additional information + % and returns it in a cell array. + fid = fopen([filename '.txt'], 'w'); + + pr = m_i.poisson; + ym = m_i.youngs * 100000 / 1e9; + rs = m_i.regional_stress; + fric = m_i.friction; + dimensions = m_i.dimensions; + res = dimensions(3); + param_id = m_i.param_id; + + switch keyword + case 'ISIP' + keyword = 'dP'; + param = (param_id * 1e-6) - 36.824; + param2 = ''; + unit = 'MPa'; + case 'FlowRate' + param = param_id; + unit = 'm^3/d'; + unit2 = 'm^3/s'; + param2 = [' (' num2str(param/86400, 5) ' ' unit2 ')']; + case 'PumpTime' + param = param_id; + unit = 'mins'; + param2 = ''; + end + + fprintf(fid, 'Calculation for: \n'); + fprintf(fid, '%s %u %s %s \n%u m resolution \n', keyword, param, unit, param2); + fprintf(fid, 'Monte Carlo Iteration %u\n \n', mc); + fprintf(fid, 'Geological Conditions: \n'); + fprintf(fid, 'Poisson ratio %2.2f\n', pr); + fprintf(fid, 'Young''s modulus %2.2f GPa\n', ym); + fprintf(fid, 'Friction %2.2f\n\n', fric); + fprintf(fid, 'Regional Stress: %s \n\n', mat2str(rs)); + fprintf(fid, 'Lateral Respect Distance (lrd) data %2.3f MPa threshold\n', latdistthreshold); + fprintf(fid, 'Mean lrd %4.2f \n', latdist(1)*1000); + fprintf(fid, 'Median lrd %4.2f \n', latdist(2)*1000); + fprintf(fid, 'Maximum lrd %4.2f \n', latdist(3)*1000); + fprintf(fid, 'Minimum lrd %4.2f \n', latdist(4)*1000); + + fclose(fid); +end + +function [stress_max, xyz] = max_stress(x, y, z, stress) + % MAX_STRESS compute the max or min stress over all depths + + % [MS, XYZ] = MAX_STRESS(X, Y, Z, STRESS) + % Takes the X, Y, Z and STRESS vectors, representing a grid of + % points and stress and returns the maximum and minimum stresses at + % each radius in MPa + + stress = stress/10; % divide by 10 to get MPa + + % reshape the data to be a cube + stress = reshape(stress, length(unique(x)), length(unique(y)), ... + length(unique(z)), []); + stress_max = stress; + + xyz.x = x(:); + xyz.y = y(:); + xyz.z = z(:); + + if ndims(stress_max) > 2 %#ok % check how many levels of z we have + + % calculate using abs max at each x,y and make sign right + [stress_max, max_position] = max(abs(stress),[],3); + [~, min_position] = min(stress,[],3); + stress_sign = (min_position == max_position) .* -1.0 ; + stress_sign(stress_sign == 0) = 1; + stress_max = stress_max .* stress_sign; + xyz.z = z(max_position); + + end +end + +function fig = plot_stress(fracX, fracY, stress, xyz, plot_fracs, stresstype, resolution) + %PLOT_STRESS. Plots the stress maps + + % PLOT_STRESS(FRACX, FRACY, STRESS, XYZ, PLOT_FRACS) + % Reads in the maximum STRESS and XYZ and plots the stress maps, with + % the option to overlay the fractures + % + % FRACX and FRACY are the x and y coordinates of the fractures + % respectively. + % STRESS is the stress in a m x n matrix + % XYZ is a structure containing the x, y and z values as vectors + % PLOT_FRACS is logical and states whether the fractures should be + % plotted or not. + + if ~exist('plot_fracs', 'var') || isempty(plot_fracs) + plot_fracs = 1; + end + + % set the logical for plotting the fractures and the fontsize + PLOT_FRACS = logical(plot_fracs); + FONTSIZE = 14; + + % set the boundaries for clipping the data + if resolution == 50 + clip = [-4e-4 4e-4]; + interval = 1e-5; + else + clip = [-4e-5 4e-5]; + interval = 1e-6; + end + + % clip the data + stress(stress < clip(1)) = clip(1); + stress(stress > clip(2)) = clip(2); + + % set the colour levels + c_levels = clip(1):interval:clip(2); + + % generate a new figure + scrsz = get(groot,'ScreenSize'); + fig = figure('Position',[10 scrsz(4)/2-100 scrsz(3)/2 scrsz(4)/2]); + + % set the colourmap + colormap(anatolia); + + % contour plot the stress and set the fontsize of the axis + contourf(unique(xyz.x), unique(xyz.y), stress, c_levels, 'LineColor', 'none'); + shading interp + + set(gca, 'FontSize', FONTSIZE);%ax.FontSize = FONTSIZE; + + % plot contour lines - black for positive and grey for negative + hold on + contour(unique(xyz.x), unique(xyz.y), stress, 0:1e-5:1e-4, 'k'); + contour(unique(xyz.x), unique(xyz.y), stress, 0:-1e-5:-1e-4, 'LineColor', [.5 .5 .5]); + hold off + + % set the labels + xlabel('X (m)','FontSize',FONTSIZE) + ylabel('Y (m)','FontSize',FONTSIZE) + caxis (clip) + + % set axis details + axis equal + axis([min(xyz.x) max(xyz.x) min(xyz.y) max(xyz.y)]); + + + if PLOT_FRACS + fcolour = [0 .7 0]; % green + + % get the number of fractures + m = length(fracX(1,:)); + + % cycle round and plot + for k = 1:m + line(fracX(:,k), fracY(:,k), 'LineStyle', '-', 'Color', fcolour); + end + end + + % add the colourbar + if resolution == 50 + hcb = colorbar('Location', 'EastOutside',... + 'Ticks', [-4e-4,-2e-4,0,2e-4,4e-4],... + 'TickLabels', {'-0.0004', '-0.0002', '0' '0.0002', '0.0004'}, ... + 'FontSize', 14); + else + hcb = colorbar('Location', 'EastOutside',... + 'Ticks', [-4e-5,-2e-5,0,2e-5,4e-5],... + 'TickLabels', {'-0.00004', '-0.00002', '0' '0.00002', '0.00004'}, ... + 'FontSize', 14); + end + title(hcb, [stresstype ' Stress (MPa)']); + %set(get(get(hcb, 'Label'), 'String'), [stresstype ' Stress (MPa)']); + + + % set the renderer + set(gcf,'renderer','painters') + +end + +function [a,b,r,s] = conv_hull(stress, x, y) + %CONV_HULL. Obtains all the points of the max stress v radius with the conve hull + %bound + % + % [A, B, R, S] = CONV_HULL(STRESS, X, Y) + % + % STRESS is the stress in a mxn matrix + % X is the x values as vectors + % Y is the y values as vectors + % + % A and B are the parameters to define the convex hull + % R is the radius in km and S is the max stress + + + % set the centre point + pt = [0,0]; + + % calculate the radius from the centre of the fractures and sort + radius = sqrt((x-pt(1)).^2 + (y-pt(2)).^2); + [radius, idx]= sort(radius); + + % sort the stress into the same order as the radius + stress = stress(:); + stress = stress(idx); + + % get the unique radius + r = unique(radius); + + % create a new array for the max stress at each radius + s = zeros(size(r)); + + % loop round to populate s + for n = 1:length(s) + ind = radius == r(n); + s(n) = max(stress(ind)); + end + + % find all values of stress less than 0 and get rid of them + ind = s <= 0; + r(ind) = []; + s(ind) = []; + + % Get rid of any radius close in (<0.25 km) + r1 = r(r>0.25); + s1 = s(r>0.25); + + % generate the convexhull + ind = convhull(log(r1), log(s1)); + + % number of points + i = numel(ind); + + % loop round to find the first point on the curve + while (r1(ind(i))<0.750) %& i>1 + i = i-1; + end + + R1 = r1(ind(i)); + E1 = s1(ind(i)); + R2 = r1(ind(i+1)); + E2 = s1(ind(i+1)); + b = log(E1/E2)/log(R1/R2); + a = E1/R1^b; +end + + +function descr = outputText(m_i, keyword, latdist, latdistthreshold, mc) + % Function which generates additional information + % and returns it in a cell array. + pr = m_i.poisson; + ym = m_i.youngs * 100000 / 1e9; + rs = m_i.regional_stress; + fric = m_i.friction; + dimensions = m_i.dimensions; + res = dimensions(3); + param_id = m_i.param_id; + + switch keyword + case 'ISIP' + keyword = 'dP'; + param = (param_id * 1e-6) - 36.824; + param2 = ''; + unit = 'MPa'; + case 'FlowRate' + param = param_id; + unit = 'm^3/d'; + unit2 = 'm^3/s'; + param2 = [' (' num2str(param/86400, 5) ' ' unit2 ')']; + case 'PumpTime' + param = param_id; + unit = 'mins'; + param2 = ''; + end + + descr = {'Calculation for:'; + [keyword ': ' num2str(param, 5) ' ' unit param2 ', ' ... + num2str(res) 'm resolution']; + ['Monte Carlo iteration: ' num2str(mc)]; + ' '; + 'Geological Conditions:'; + ['Poisson ratio ' num2str(pr)]; + ['Young''s modulus ' num2str(ym) ' GPa']; + ['Friction ' num2str(fric)] + ''; + 'Regional Stress:'; + mat2str(rs) + ''; + ['Lateral Respect Distance (lrd) data ' num2str(latdistthreshold) ... + ' MPa threshold: '] + ['Mean lrd ' num2str(latdist(1)*1000) ' m'] + ['Median lrd ' num2str(latdist(2)*1000) ' m'] + ['Maximum lrd ' num2str(latdist(3)*1000) ' m'] + ['Minimum lrd ' num2str(latdist(1)*1000) ' m']; + }; + + +end + +function [mean_ld, median_ld, max_ld, min_ld] = lateralDistance(x, y, stress, latdistthreshold) + %LATERALDISTANCE Calculates the lateral distance using a convex + % hull + % + % [MEAN_LD, MEDIAN_LD, MAX_LD, MIN_LD] = LATERALDISTANCE(X, Y, Z, STRESS, FRAC_CENT) + % + % X is the x values as a vector + % Y is the y values as a vector + % STRESS is the stress in a m x n matrix + % LATDISTTHRESHOLD is the stress threshold to calculate the lateral + % respect distance for + + % check if we are in km, if not, convert + if max(abs(x)) > 10 + x = x/1000; + y = y/1000; + end + + % create a radius vector + radius = zeros(1,50); + + % loop round all monte carlo iterations + for jj = 1:50 + [a, b] = conv_hull(stress(:,jj), x, y); + + % compute the radius + radius(jj) = (latdistthreshold/a)^(1/b); + end + + mean_ld = mean(radius); + median_ld = median(radius); + max_ld = max(radius); + min_ld = min(radius); +end + +function cmap = anatolia() + % the Anatolia colormap + + cmap = ... + [0.21176470816135406 0.15294118225574493 0.99215686321258545; + 0.21398124098777771 0.1877237856388092 0.99232739210128784; + 0.21619778871536255 0.22250640392303467 0.99249786138534546; + 0.21841432154178619 0.25728902220726013 0.99266839027404785; + 0.22063086926937103 0.29207161068916321 0.99283885955810547; + 0.22284740209579468 0.32685422897338867 0.99300938844680786; + 0.22506394982337952 0.36163684725761414 0.99317985773086548; + 0.22728048264980316 0.39641943573951721 0.99335038661956787; + 0.22949701547622681 0.43120205402374268 0.99352091550827026; + 0.23171356320381165 0.46598467230796814 0.99369138479232788; + 0.23393009603023529 0.5007672905921936 0.99386191368103027; + 0.23614664375782013 0.53554987907409668 0.99403238296508789; + 0.23836317658424377 0.57033246755599976 0.99420291185379028; + 0.24057972431182861 0.60511511564254761 0.9943733811378479; + 0.24279625713825226 0.63989770412445068 0.99454391002655029; + 0.2450128048658371 0.67468029260635376 0.99471437931060791; + 0.24722933769226074 0.70946294069290161 0.9948849081993103; + 0.24944587051868439 0.74424552917480469 0.9950554370880127; + 0.25166240334510803 0.77902811765670776 0.99522590637207031; + 0.25387895107269287 0.81381076574325562 0.99539643526077271; + 0.25609549880027771 0.84859335422515869 0.99556690454483032; + 0.25831204652786255 0.88337594270706177 0.99573743343353271; + 0.260528564453125 0.91815859079360962 0.99590790271759033; + 0.26274511218070984 0.9529411792755127 0.99607843160629272; + 0.35490196943283081 0.958823561668396 0.99656862020492554; + 0.44705882668495178 0.96470588445663452 0.99705880880355835; + 0.53921568393707275 0.970588207244873 0.99754899740219116; + 0.63137257099151611 0.97647058963775635 0.99803924560546875; + 0.7235293984413147 0.98235297203063965 0.99852943420410156; + 0.81568628549575806 0.98823529481887817 0.99901962280273438; + 0.90784311294555664 0.9941176176071167 0.99950981140136719; + 1 1 1; + 1 1 1; + 0.99901962280273438 0.99215686321258545 0.91372549533843994; + 0.99803924560546875 0.9843137264251709 0.82745099067687988; + 0.99705880880355835 0.97647058963775635 0.74117648601531982; + 0.99607843160629272 0.9686274528503418 0.65490198135375977; + 0.9950980544090271 0.96078431606292725 0.56862747669219971; + 0.9941176176071167 0.9529411792755127 0.48235294222831726; + 0.99313724040985107 0.94509804248809814 0.3960784375667572; + 0.99215686321258545 0.93725490570068359 0.30980393290519714; + 0.989599347114563 0.89650470018386841 0.29991474747657776; + 0.98704177141189575 0.85575449466705322 0.29002559185028076; + 0.98448425531387329 0.815004289150238 0.28013640642166138; + 0.981926679611206 0.77425402402877808 0.27024725079536438; + 0.97936916351318359 0.73350381851196289 0.260358065366745; + 0.97681158781051636 0.69275361299514771 0.25046887993812561; + 0.9742540717124939 0.65200340747833252 0.24057972431182861; + 0.97169649600982666 0.61125320196151733 0.23069053888320923; + 0.9691389799118042 0.57050299644470215 0.22080136835575104; + 0.966581404209137 0.529752790927887 0.21091219782829285; + 0.9640238881111145 0.48900255560874939 0.20102302730083466; + 0.96146631240844727 0.4482523500919342 0.19113385677337646; + 0.9589087963104248 0.407502144575119 0.18124467134475708; + 0.95635122060775757 0.36675190925598145 0.17135550081729889; + 0.95379370450973511 0.32600170373916626 0.1614663302898407; + 0.95123612880706787 0.28525149822235107 0.15157715976238251; + 0.94867861270904541 0.24450127780437469 0.14168798923492432; + 0.94612103700637817 0.20375107228755951 0.13179880380630493; + 0.94356352090835571 0.16300085186958313 0.12190964072942734; + 0.94100594520568848 0.12225063890218735 0.11202046275138855; + 0.938448429107666 0.081500425934791565 0.10213129222393036; + 0.93589085340499878 0.040750212967395782 0.092242114245891571; + 0.93333333730697632 0 0.08235294371843338]; +end \ No newline at end of file