diff --git a/autodiff/test-suite/TestCase.m b/autodiff/test-suite/TestCase.m index f6ab7dddb..339bfa29e 100644 --- a/autodiff/test-suite/TestCase.m +++ b/autodiff/test-suite/TestCase.m @@ -101,8 +101,8 @@ test.toolbarOptions([1:2:end, 2:2:end]) = horzcat(names, values); % Set plotting function if not given if isempty(test.plotFn) - test.plotFn = @(G, v, varargin) ... - plotToolbar(G, v, test.toolbarOptions{:}, varargin{:}); + test.plotFn = @(G, v, cells, varargin) ... + plotToolbar(G, v, cells, test.toolbarOptions{:}, varargin{:}); end % Compute test case ID if opt.computeID @@ -278,12 +278,20 @@ require mrst-gui if nargin == 1 || ischar(varargin{1}) - v = test.model.rock; + if isa(test.model, 'CompositeModel') + assert(isfield(test.model.submodels, 'Reservoir'), ... + ['Did not find submodel called ''Reservoir'', ', ... + 'please provide values for plotting.']); + v = test.model.submodels.Reservoir.rock; + else + v = test.model.rock; + end else v = varargin{1}; varargin = varargin(2:end); end opt = struct('Name' , '' , ... + 'cells' , [] , ... 'plotWells', true, ... 'wellOpts' , {{}}, ... 'camlight' , true); @@ -294,7 +302,10 @@ end h = test.figure('Name', Name); G = test.getVisualizationGrid(); - test.plotFn(G, v, varargin{:}); + if isempty(opt.cells) + opt.cells = true(G.cells.num, 1); + end + test.plotFn(G, v, opt.cells, varargin{:}); if opt.plotWells test.plotWells(opt.wellOpts{:}); end @@ -322,47 +333,55 @@ function varargout = plotWells(test, varargin) % Plot wells in the model - if ~isfield(test.schedule.control(1), 'W') + W = []; + if isa(test.model, 'CompositeModel') + if isfield(test.model.submodels, 'Wellbore') + W = test.model.submodels.Wellbore.wells; + else + warning(['Did not find submodel called ''Wellbore'', ', ... + 'cannot plot wells.']); + end + elseif isfield(test.schedule.control(1), 'W') + W = test.schedule.control(1).W; + end + if isempty(W) return; end - W = test.schedule.control(1).W; - if ~isempty(W) - G = test.getVisualizationGrid(); - if G.griddim == 3 - ax = gca; dz = ax.ZLim(2) - ax.ZLim(1); - varargout = cell(nargout, 1); - [varargout{:}] = plotWell(G, W, 'color' , 'k' , ... - 'height', 0.15*dz, ... - varargin{:} ); - else - opt = struct('fontsize', 12 , ... - 'color' , 'k', ... - 'color2' , 'k'); - [opt, varargin] = merge_options(opt, varargin{:}); - nw = numel(W); - [hw, htext] = deal(zeros(nw,1)); - gray = [1,1,1]*0.8; - hold on - for w = 1:numel(W) - color = opt.color; - if W(w).sign < 0, color = opt.color2; end - x = G.cells.centroids(vertcat(W(w).cells), :); - hw(w) = plot(x(1), x(2) , 'o' , ... - 'Color' , color, ... - 'markerFaceColor', gray , ... - 'markerSize' , 8 , ... - varargin{:} ); - if opt.fontsize > 0 - htext(w) = text(x(1), x(2), W(w).name, ... - 'FontSize' , opt.fontsize, ... - 'VerticalAlignment', 'bottom' , ... - 'interp' , 'latex' , ... - 'Color' , color ); - end + G = test.getVisualizationGrid(); + if G.griddim == 3 + ax = gca; dz = ax.ZLim(2) - ax.ZLim(1); + varargout = cell(nargout, 1); + [varargout{:}] = plotWell(G, W, 'color' , 'k' , ... + 'height', 0.15*dz, ... + varargin{:} ); + else + opt = struct('fontsize', 12 , ... + 'color' , 'k', ... + 'color2' , 'k'); + [opt, varargin] = merge_options(opt, varargin{:}); + nw = numel(W); + [hw, htext] = deal(zeros(nw,1)); + gray = [1,1,1]*0.8; + hold on + for w = 1:numel(W) + color = opt.color; + if W(w).sign < 0, color = opt.color2; end + x = G.cells.centroids(vertcat(W(w).cells), :); + hw(w) = plot(x(1), x(2) , 'o' , ... + 'Color' , color, ... + 'markerFaceColor', gray , ... + 'markerSize' , 8 , ... + varargin{:} ); + if opt.fontsize > 0 + htext(w) = text(x(1), x(2), W(w).name, ... + 'FontSize' , opt.fontsize, ... + 'VerticalAlignment', 'bottom' , ... + 'interp' , 'latex' , ... + 'Color' , color ); end - hold off - varargout = {hw, htext}; end + hold off + varargout = {hw, htext}; end end @@ -371,8 +390,20 @@ function ok = hasDrivingForce(test, force) % Check if the test case has a given driving force - ok = isfield(test.schedule.control(1), force) ... - && ~isempty(test.schedule.control(1).(force)); + if isa(test.model, 'CompositeModel') + switch force + case 'W' + mdl = 'Wellbore'; + case {'bc', 'src'} + mdl = 'Reservoir'; + end + ok = isfield(test.model.submodels, mdl) ... + && isfield(test.schedule.control(1).(mdl), force) ... + && ~isempty(test.schedule.control(1).(mdl).(force)); + else + ok = isfield(test.schedule.control(1), force) ... + && ~isempty(test.schedule.control(1).(force)); + end end diff --git a/autodiff/test-suite/setup-functions/layered_nls_test_wog.m b/autodiff/test-suite/setup-functions/layered_nls_test_wog.m deleted file mode 100644 index 314d73bf4..000000000 --- a/autodiff/test-suite/setup-functions/layered_nls_test_wog.m +++ /dev/null @@ -1 +0,0 @@ -% diff --git a/modules/geothermal/examples/geothermalExampleWellbore.m b/modules/geothermal/examples/geothermalExampleWellbore.m index 2ee208e77..70eb6fd2d 100644 --- a/modules/geothermal/examples/geothermalExampleWellbore.m +++ b/modules/geothermal/examples/geothermalExampleWellbore.m @@ -52,10 +52,7 @@ test.schedule.control(2).Wellbore = ctrl; %% Set up linear solver -lsol = setUpGeothermalWellboreLinearSolver(test.model); -test.model.G.cells.num = test.model.submodels.Reservoir.G.cells.num + ... -nnz(test.model.submodels.Wellbore.G.cells.type == 0 & ... - test.model.submodels.Wellbore.G.cells.hybrid == 0); +[lsol, test.model] = setUpGeothermalWellboreLinearSolver(test.model); %% Simulate problem = test.getPackedSimulationProblem('LinearSolver', lsol); diff --git a/modules/geothermal/models/wellbore/WellboreModel.m b/modules/geothermal/models/wellbore/WellboreModel.m index 79f7e1838..78ee93071 100644 --- a/modules/geothermal/models/wellbore/WellboreModel.m +++ b/modules/geothermal/models/wellbore/WellboreModel.m @@ -489,7 +489,6 @@ function [eqs, names, types, state] = getMassFluxEquations(model, state) % Equations realting segment mass fluxes to the pressure by means % of a wellbore friction model - dp = model.getProp(state, 'FrictionLoss'); [p, s, rho] = model.parentModel.getProps(state, ... @@ -497,6 +496,7 @@ 's' , ... 'Density' ... ); + v = model.getProp(state, 'massFlux'); nph = model.getNumberOfPhases(); @@ -511,10 +511,20 @@ g = norm(model.parentModel.gravity); dz = model.parentModel.operators.Grad(model.G.cells.centroids(:,3)); dpw = model.parentModel.operators.Grad(p); + pot = dpw - rhoMix.*g.*dz; + + % Assemble equation. We add the mass flux so that this reduces to a + % Darcy-type equation for very low friction loss to avoid singular + % systems. + % Convert mass flux to velocity + [Di, Do] = deal(0, model.G.faces.radius.*2); + if size(Do, 2) == 2 + Di = Do(:,1); + Do = Do(:,2); + end + v = v./(pi*rhoMix.*((Do/2).^2 - (Di/2).^2)); + eqs = (pot - dp) + 10*v; - pot = dpw - rhoMix.*g.*dz; - - eqs = (pot - dp)./(1*atm); eqs = {eqs}; names = {'flux'}; types = {'face'}; @@ -676,7 +686,7 @@ % re-initialize the values and change the controls so that the next % step keeps within the prescribed ranges. - + if ~isfield(state, 'wells'), state.wells = model.wells; end if ~isfield(state, 'groups'), state.groups = model.groups; end % Get driving forces and check that it's non-empty diff --git a/modules/geothermal/utils/addHybridCellsToWell.m b/modules/geothermal/utils/addHybridCellsToWell.m index 82863a83a..efcb535e9 100644 --- a/modules/geothermal/utils/addHybridCellsToWell.m +++ b/modules/geothermal/utils/addHybridCellsToWell.m @@ -3,11 +3,10 @@ % % SYNOPSIS: % W = addHybridCellsToWell(W, G, rock) -% setup = spe10_wo(W, G, rock, 'pn1', pv1, ...) % % DESCRIPTION: -% For each well in W, this function finds all hybrid cells adjacent ot -% the cells of the well and adds them to the well. Hybrid are grid faces +% For each well in W, this function finds all hybrid cells adjacent to the +% cells of the well and adds them to the well. Hybrid cells are grid faces % that are assigned a volume in the computational grid. This function is % useful if the wells are added using a function that does not recognize % hybrid cells, e.g., addWellFromTrajectory. @@ -79,6 +78,7 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). function W = addWellCells(W, G, rock, cells, aperture, varargin) % Add cells + cells0 = W.cells; nc0 = numel(W.cells); W.cells = [W.cells; cells]; nc = numel(W.cells); @@ -103,6 +103,16 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). W.cells = W.cells(ix); W.WI = W.WI(ix); W.dZ = W.dZ(ix); + + if isfield(W, 'trajectory') + assert(all(W.trajectory.cell == cells0), ... + 'W.trajectory.cell must match W.cells'); + W.trajectory.cell = W.cells; + vec = [W.trajectory.vec; [zeros(numel(cells), 2), dz]]; + W.trajectory.vec = vec(ix,:); + weight = [W.trajectory.weight; ones(numel(cells), 1)]; + W.trajectory.weight = weight(ix); + end % Add remaining fields by repeating W.(fn)(1) for fn = fieldnames(W)' diff --git a/modules/geothermal/utils/convertToReservoirWellboreModel.m b/modules/geothermal/utils/convertToReservoirWellboreModel.m index cd3588ccc..5166b169e 100644 --- a/modules/geothermal/utils/convertToReservoirWellboreModel.m +++ b/modules/geothermal/utils/convertToReservoirWellboreModel.m @@ -37,11 +37,15 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). % Optional input arguments opt = struct('name', []); [opt, varargin] = merge_options(opt, varargin{:}); - if isempty(opt.name), opt.name = [setup.name, '_wb']; end + if isfield(setup, 'name') && isempty(opt.name), opt.name = [setup.name, '_wb']; end % Get reservoir model and construct wellbore model reservoirModel = setup.model; - W = setup.schedule.control(1).W; + if isfield(setup, 'W') + W = setup.W; + else + W = setup.schedule.control(1).W; + end wellboreModel = WellboreModel(reservoirModel, W, varargin{:}); % Make composite model from the two model = CompositeModel({reservoirModel, wellboreModel}, ... @@ -70,36 +74,38 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). state0.Wellbore.s = state0.Wellbore.s(wc); state0.Wellbore.T = state0.Wellbore.T(wc); - % Convert schedule - schedule = setup.schedule; - control = schedule.control; - for i = 1:numel(control) - [rctrl, wctrl] = deal(control(i)); - [rctrl.W, rctrl.groups] = deal([]); - if ~isfield(rctrl, 'bc' ), rctrl.bc = []; end - if ~isfield(rctrl, 'src'), rctrl.src = []; end - [wctrl.bc, wctrl.src] = deal([]); - schedule.control(i).Reservoir = rctrl; - schedule.control(i).Wellbore = wctrl; - end - if isfield(schedule.control, 'W') - schedule.control = rmfield(schedule.control, 'W'); - end - if isfield(schedule.control, 'src') - schedule.control = rmfield(schedule.control, 'src'); - end - if isfield(schedule.control, 'bc') - schedule.control = rmfield(schedule.control, 'bc'); - end - if isfield(schedule.control, 'groups') - schedule.control = rmfield(schedule.control, 'groups'); - end - % Update setup setup.state0 = state0; setup.model = model; - setup.schedule = schedule; - setup.name = opt.name; + if ~isempty(opt.name), setup.name = opt.name; end + + if isfield(setup, 'schedule') || isa(setup, 'TestCase') + % Convert schedule + schedule = setup.schedule; + control = schedule.control; + for i = 1:numel(control) + [rctrl, wctrl] = deal(control(i)); + [rctrl.W, rctrl.groups] = deal([]); + if ~isfield(rctrl, 'bc' ), rctrl.bc = []; end + if ~isfield(rctrl, 'src'), rctrl.src = []; end + [wctrl.bc, wctrl.src] = deal([]); + schedule.control(i).Reservoir = rctrl; + schedule.control(i).Wellbore = wctrl; + end + if isfield(schedule.control, 'W') + schedule.control = rmfield(schedule.control, 'W'); + end + if isfield(schedule.control, 'src') + schedule.control = rmfield(schedule.control, 'src'); + end + if isfield(schedule.control, 'bc') + schedule.control = rmfield(schedule.control, 'bc'); + end + if isfield(schedule.control, 'groups') + schedule.control = rmfield(schedule.control, 'groups'); + end + setup.schedule = schedule; + end % Set visualization grid setup.visualizationGrid = model.submodels.Reservoir.G; diff --git a/modules/geothermal/utils/initializeGeothermalEquilibrium.m b/modules/geothermal/utils/initializeGeothermalEquilibrium.m index 09219763e..20b827426 100644 --- a/modules/geothermal/utils/initializeGeothermalEquilibrium.m +++ b/modules/geothermal/utils/initializeGeothermalEquilibrium.m @@ -60,8 +60,8 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). % Set up and solve ODE defining geothermal equilibrilum %---------------------------------------------------------------------% - [dpdz, dTdz, p0, T0, z, z0] = getPhysicalQuantities(model, opt); - [p,T,sol] = solveODE(dpdz, dTdz, p0, T0, z, z0); + [dpdz, dTdz, p0, T0, zc, zb, z0] = getPhysicalQuantities(model, opt); + [p,T,sol] = solveODE(dpdz, dTdz, p0, T0, zc, zb, z0); %---------------------------------------------------------------------% end @@ -86,7 +86,7 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). end %-------------------------------------------------------------------------% -function [dpdz, dTdz, p0, T0, z, z0] = getPhysicalQuantities(model, opt) +function [dpdz, dTdz, p0, T0, zc, zb, z0] = getPhysicalQuantities(model, opt) % Get physical quantities defining the ODE % Pressure gradient due to density @@ -99,13 +99,15 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). p0 = opt.datumPressure; T0 = opt.datumTemperature; % Verical coordinates in the grid - z = model.G.cells.centroids(:,3); + zc = model.G.cells.centroids(:,3); + bfaces = boundaryFaces(model.G); + zb = model.G.faces.centroids(bfaces,3); z0 = opt.datumDepth; end %-------------------------------------------------------------------------% -function [p,T,sol] = solveODE(dpdz, dTdz, p0, T0, z, z0) +function [p,T,sol] = solveODE(dpdz, dTdz, p0, T0, zc, zb, z0) % Solve ODE to find equilibrium pressure and temperature % Set initial condition @@ -116,6 +118,7 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). odeopts = odeset('AbsTol', 1.0e-10, 'RelTol', 5.0e-8); % Solve ODE [solUp, solDown] = deal([]); + z = [zc; zb]; if z0 > min(z) zlimUp = [z0, min(z)]; solUp = ode45(dydz, zlimUp, y0, odeopts); @@ -126,7 +129,7 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). end assert(~(isempty(solUp) && isempty(solDown)), 'Something went wrong'); % Evaluate in cells - [p,T] = evaluateFn(solUp, solDown, z0, z); + [p,T] = evaluateFn(solUp, solDown, z0, zc); sol = @(z) evaluateFn(solUp, solDown, z0, z); end diff --git a/modules/geothermal/utils/processWellboreTrajectories.m b/modules/geothermal/utils/processWellboreTrajectories.m index 0049554eb..bfedfec46 100644 --- a/modules/geothermal/utils/processWellboreTrajectories.m +++ b/modules/geothermal/utils/processWellboreTrajectories.m @@ -71,11 +71,15 @@ This file is part of The MATLAB Reservoir Simulation Toolbox (MRST). trajectories = cell(nw,1); for i = 1:nw if wellsGiven - W(i) = topoSortWellCells(G, W(i)); %#ok - coords{i} = G.cells.centroids(W(i).cells,:); - traj = computeTraversedCells(G, coords{i}); - if isfield(G.cells, 'hybrid') && any(G.cells.hybrid(W(i).cells)) - traj = addHybridCells(G, W(i), traj); + if isfield(W(i), 'trajectory') + traj = W(i).trajectory; + else + W(i) = topoSortWellCells(G, W(i)); %#ok + coords{i} = G.cells.centroids(W(i).cells,:); + traj = computeTraversedCells(G, coords{i}); + if isfield(G.cells, 'hybrid') && any(G.cells.hybrid(W(i).cells)) + traj = addHybridCells(G, W(i), traj); + end end trajectories{i} = traj; else diff --git a/modules/geothermal/utils/setUpGeothermalWellboreLinearSolver.m b/modules/geothermal/utils/setUpGeothermalWellboreLinearSolver.m index 1bfbd2858..21126aac9 100644 --- a/modules/geothermal/utils/setUpGeothermalWellboreLinearSolver.m +++ b/modules/geothermal/utils/setUpGeothermalWellboreLinearSolver.m @@ -1,8 +1,12 @@ -function solver = setUpGeothermalWellboreLinearSolver(model, varargin) +function [solver, model] = setUpGeothermalWellboreLinearSolver(model, varargin) %Set up a linear solver for composite geothermal models with WellboreModel % Optional input arguments - opt = struct('useCPR', true); + opt = struct( ... + 'useCPR', true, ... + 'tolerance', 1e-4, ... + 'maxIterations', 200 ... + ); opt = merge_options(opt, varargin{:}); if ~opt.useCPR @@ -15,8 +19,8 @@ solver.decoupling = 'quasiimpes'; end % Set maximum iterations and tolerance - solver.maxIterations = 200; - solver.tolerance = 1e-4; + solver.maxIterations = opt.maxIterations; + solver.tolerance = opt.tolerance; % Set block size bz = 1 + model.submodels.Reservoir.thermal; @@ -41,7 +45,9 @@ solver.keepNumber = bz*(ncr + ncw); solver.amgcl_setup.cell_size = ncr + ncw; solver.reduceToCell = false; - + + model.G.cells.num = ncr + ncw; + end %{ diff --git a/modules/upr/pebi2D/lineSites2D.m b/modules/upr/pebi2D/lineSites2D.m index f9fecd46e..701cf719a 100644 --- a/modules/upr/pebi2D/lineSites2D.m +++ b/modules/upr/pebi2D/lineSites2D.m @@ -143,6 +143,7 @@ end wellSpace = sqrt(sum(diff(p,1,1).^2,2)); wellSpace = min(wellSpace([1 1:end]), wellSpace([1:end end])); + if size(wellSpace, 2) > 1, wellSpace = reshape(wellSpace, [], 1); end end keep = 1:size(p,1);