Errors Resulting from Geothermal and EDFM Model Operators Setup

[David E]

Hello there,

I will appreciate some help to resolve errors arising from combining EDFM model with geothermal module.

I have a script that runs an EDFM model and initially simulates for pressure in fractures. The second part of the script contains the geothermal component to account for both pressure and heat flow in the same EDFM model.

The first issue is that the model operators setup (setupEDFMOperatorsTPFA) for only pressure simulation results in the following error when evaluating molecular transmissibility in the geothermal model:

% setup model

     model = GeothermalModel(G, G.rock, fluid);

      model.operators = setupEDFMOperatorsTPFA(G, G.rock, tol);

 Solving timestep 1/1:               -> 1 Millisecond

The logical indices contain a true value outside of the array bounds.

Error in ComponentTotalDiffusiveFlux>@(T)T(model.operators.internalConn) (line 20)
T(ix) = cellfun(@(T) T(model.operators.internalConn), T(ix), 'UniformOutput', false);
Error in ComponentTotalDiffusiveFlux/evaluateOnDomain (line 20)
T(ix) = cellfun(@(T) T(model.operators.internalConn), T(ix), 'UniformOutput', false);
...

I resolved this error by using setupOperatorsTPFA instead and ensured that the transmissibility and neighborship of internal faces are passed to the setup as shown below. 

%The transmissibility for each face in the combined grid and each NNC is computed and stored in the vector T.

     T = computeTrans(G, G.rock);
     cf = G.cells.faces(:,1);
     nf = G.faces.num;
     T = 1 ./ accumarray(cf, 1./T, [nf, 1]);
     T = [T;G.nnc.T];
     N = getNeighbourship(G, 'topological', true);
     intx = all(N ~= 0, 2);

     model.operators = setupOperatorsTPFA(G, G.rock, 'trans', T(intx),...
     'neighbors', N(intx, :));
     model.operators.T_all = T;

This work around normally allows my simulations to run without problems. However, for this model, I had to add model.operators.T_all  line for the simulation to run. Despite running it to the end without any indication of additional problems, no pressure or flux is computed in the resulting states. 

Without the model.operators.T_all  line, the simulation returns the following error stating that the assigned BC faces are outside the accounted transmissibility elements:

Solving timestep 1/1:               -> 1 Millisecond

Index exceeds the number of array elements (52545).
Error in getBoundaryConditionFluxesAD (line 63)
T = model.operators.T_all(bc.face);
Error in computeSourcesAndBoundaryConditionsAD (line 82)
[qVolBC, BCTocellMap, bcCells, qResBC] =
getBoundaryConditionFluxesAD(model, pressure, s, mob, rho, b,
forces.bc);
...

I will appreciate your suggestions on how to correctly fix these errors.

Kind regards,

David Egya

[Øystein S. Klemetsdal]

Dear David,

Thank you for pointing this out.

The reason for the behavior is that permeability and thermal conductivities are allowed to depend on pressure and transport, which means that two-point transmissibilities must be recomputed at each nonlinear iteration.

This does make the module less flexible in combining with other functionality in MRST, and we have a planned fix for the next release.

Cheers,

Øystein Klemetsdal // Research Scientist, SINTEF Digital

[David E]

Thanks, Øystein, for your response. We appreciate your hard work in making MRST modules more accessible and user-friendly. 

For the case in question, I still don't get any updated pressure and flux in the simulated states even though the model is running. Any handy recommendation on the best way to work around this issue pending the planned fix for the next release?

I'm happy to share my scripts and any other information that may be required.

Thanks in advance.

David Egya (Heriot-Watt University)  

[Øystein S. Klemetsdal]

Dear David,

Thank you for the kind words!

If you send me a script reproducing the problem, I can have a look at it.

Cheers,

- Øystein