-
Notifications
You must be signed in to change notification settings - Fork 9
/
trapsJohansen.m
124 lines (106 loc) · 4.3 KB
/
trapsJohansen.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
%% Trapping on Johansen grids
% The Johansen formation has been proposed as a potential injection site
% for CO2, in particular when it was planned to capture CO2 from the
% gas-power plant at Mongstad. A simplified representation of the Johansen
% formation was also used as a benchmark case in the Stuttgart code
% comparison study based on a paper by Eigestad et al. Herein, we consider
% the same injection point as suggested by Eigestad et al. and evaluate the
% potential for structural trapping that can be estimated from two
% different models: (i) a sector model developed by the Norwegian Petroleum
% Directorate, which was used to produce the Stuttgart benchmark test, and
% (ii) a model of a somewhat larger region derived from the CO2 Storage
% Atlas for the Norwegian North Sea.
%
%%Author : Salman Karim
%
mrstModule add co2lab libgeometry opm_gridprocessing coarsegrid matlab_bgl
%% Load NPD data: sector model
try
jdir = fullfile(mrstPath('co2lab'), 'data', 'johansen');
sector = 'NPD5';
sector = fullfile(jdir, sector);
grdecl = readGRDECL([sector '.grdecl']);
catch me
disp(' -> Download data from: http://www.sintef.no/Projectweb/MatMoRA/')
disp([' Putting data in ', jdir]);
unzip('http://www.sintef.no/project/MatMoRA/Johansen/NPD5.zip', jdir);
grdecl = readGRDECL([sector '.grdecl']);
end
% Extract the part that represents the Johansen formation
grdecl = cutGrdecl(grdecl,[1 grdecl.cartDims(1); 1 grdecl.cartDims(2); 6 11]);
try
Gs = processgrid(grdecl);
Gs = mcomputeGeometry(Gs);
catch
Gs = processGRDECL(grdecl);
Gs = computeGeometry(Gs);
end
Gts = topSurfaceGrid(Gs);
% Get the position of the well (data given from 'Sector5_Well.txt');
wi = find(Gs.cells.indexMap==sub2ind(Gs.cartDims, 48, 48, 1));
wcent = Gs.cells.centroids(wi,:);
d = sqrt(sum(bsxfun(@minus, Gts.cells.centroids, wcent(1:2)).^2, 2));
[~,wi_s] = min(d);
%% Load NPD data: full-field model
try
jdir = fullfile(mrstPath('co2lab'), 'data', 'johansen');
sector = 'FULLFIELD_IMAXJMAX';
sector = fullfile(jdir, sector);
grdecl = readGRDECL([sector '.GRDECL']);
catch me
disp(' -> Download data from: http://www.sintef.no/Projectweb/MatMoRA/')
disp([' Putting data in ', jdir]);
unzip('https://www.sintef.no/project/MatMoRA/Johansen/FULLFIELD_Eclipse.zip', jdir);
grdecl = readGRDECL([sector '.GRDECL']);
end
% Extract the part that represents the Johansen formation
grdecl = cutGrdecl(grdecl,[1 grdecl.cartDims(1); 1 grdecl.cartDims(2); 10 14]);
try
Gf = processgrid(grdecl);
Gf = mcomputeGeometry(Gf);
catch
Gf = processGRDECL(grdecl);
Gf = computeGeometry(Gf);
end
Gtf = topSurfaceGrid(Gf);
% Get the position of the well
d = sqrt(sum(bsxfun(@minus, Gtf.cells.centroids, wcent(1:2)).^2, 2));
[~,wi_f] = min(d);
%% Load atlas data
grdecl = getAtlasGrid('Johansenfm');
try
Ga = processgrid(grdecl{1});
Ga = mcomputeGeometry(Ga);
catch
Ga = processGRDECL(grdecl{1});
Ga = computeGeometry(Ga);
end
Gta = topSurfaceGrid(Ga);
% Get the position of the well
d = sqrt(sum(bsxfun(@minus, Gta.cells.centroids, wcent(1:2)).^2, 2));
[~,wi_a] = min(d);
%% Plot the three data sets
clf
zm = min(Ga.nodes.coords(:,3));
zM = max(Ga.nodes.coords(:,3));
plotCellData(Ga,Ga.cells.centroids(:,3),'FaceAlpha',.95,'EdgeColor','none');
plotGrid(Gf,'FaceColor','none','EdgeAlpha',.2,'EdgeColor','r');
plotGrid(Gs,'FaceColor','none','EdgeAlpha',.4,'EdgeColor','k');
axis tight, view(-62,60);
light, lighting phong,
light('Position',[max(Gta.cells.centroids) 4*zM],'Style','infinite');
hold on; plot3(wcent([1 1],1),wcent([1 1],2),[zm zM],'b','LineWidth',2);
legend('Atlas','Full field','Sector','Inj.pt',...
'Location','SouthOutside','Orientation','horizontal');
%% Interactive trapping: atlas grid
interactiveTrapping(Gta, 'method', 'cell', 'light', true, ...
'spillregions', true, 'colorpath', false, 'injpt', wi_a);
view(-80,64);
%% Interactive trapping: NPD sector grid
interactiveTrapping(Gts, 'method', 'cell', 'light', true, ...
'spillregions', true, 'colorpath', false, 'injpt', wi_s);
view(-80,64);
%% Interactive trapping: NPD full-field grid
interactiveTrapping(Gtf, 'method', 'cell', 'light', true, ...
'spillregions', true, 'colorpath', false, 'injpt', wi_f);
view(-80,64);