Iterative multiscale gradient computation for heterogeneous subsurface flow
Jesus de Moraes, Rafael; de Zeeuw, Wessel; R. P. Rodrigues, José; Hajibeygi, Hadi; Jansen, Jan Dirk
DOI
10.1016/j.advwatres.2019.05.016
Publication date
2019
Document Version
Final published version
Published in
Advances in Water Resources
Citation (APA)
Jesus de Moraes, R., de Zeeuw, W., R. P. Rodrigues, J., Hajibeygi, H., & Jansen, J. D. (2019). Iterative
multiscale gradient computation for heterogeneous subsurface flow. Advances in Water Resources, 129,
210-221. https://doi.org/10.1016/j.advwatres.2019.05.016
Important note
To cite this publication, please use the final published version (if applicable).
Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
This work is downloaded from Delft University of Technology.
ContentslistsavailableatScienceDirect
Advances
in
Water
Resources
journalhomepage:www.elsevier.com/locate/advwatres
Iterative
multiscale
gradient
computation
for
heterogeneous
subsurface
flow
Rafael
J.
de
Moraes
a,b,∗,
Wessel
de
Zeeuw
a,
José Roberto
P.
Rodrigues
b,
Hadi
Hajibeygi
a,
Jan
Dirk
Jansen
aa Department of Geoscience and Engineering, Faculty of Civil Engineering and Geosciences, Delft University of Technology, P.O. Box 5048, Delft 2600, the Netherlands b Petrobras Research and Development Center CENPES, Av. Horácio Macedo 950, Cidade Universitária, Rio de Janeiro, RJ 21941-915, Brazil
a
r
t
i
c
l
e
i
n
f
o
Keywords:
Subsurface flow Gradient computation Multiscale methods
Iterative multiscale finite volume Direct method
Adjoint method
a
b
s
t
r
a
c
t
Weintroduceasemi-analyticaliterativemultiscalederivativecomputationmethodologythatallowsforerror controlandreductiontoanydesiredaccuracy,uptofine-scaleprecision.Themodelresponsesarecomputedby themultiscaleforwardsimulationofflowinheterogeneousporousmedia.Thederivativecomputationmethod isbasedontheaugmentationofthemodelequationandstatevectorswiththesmoothingstagedefinedbythe iterativemultiscalemethod.Intheformulation,weavoidadditionalcomplexityinvolvedincomputingpartial derivativesassociatedtothesmoothingstep.Weaccountforitasanapproximatederivativecomputationstage. Thenumericalexperimentsillustratehowthenewlyintroducedderivativemethodcomputesmisfitobjective functiongradientsthatconvergetofine-scaleoneastheiterativemultiscaleresidualconverges.Therobustness ofthemethodologyisinvestigatedfortestcaseswithhighcontrastpermeabilityfields.Theiterativemultiscale gradientmethodcastsapromisingapproach,withminimalaccuracy-efficiencytradeoff,forlarge-scale hetero-geneousporousmediaoptimizationproblems.
1. Introduction
Derivativecomputationisanimportantaspectofgradient-based op-timizationalgorithms.Whentheobjective(orcost)functionevaluation involvesthenumericalsimulationofdiscretizedpartialderivative equa-tions(PDE),efficientgradientcomputationisofutmostimportance.It iswelldocumentedintheliteraturethatthemostefficientand accu-ratewayofcomputingderivativesaretheanalyticaldirect(Anterion etal.,1989;Rodrigues,2006;Oliveretal.,2008)(whenthenumberof costfunctionalsisgreaterthanthenumberofparameters)andAdjoint (Chaventetal.,1975;Lietal.,2003;Oliveretal.,2008;Kraaijevanger etal.,2007;Rodrigues,2006;Jansen,2011)(ifthenumberof param-etersisgreaterthanthenumbercostfunctionals)methods.However, evenwhenconsideringefficient gradientmethods, duetothe neces-sitytoevaluatetheforwardmodelandthederivativeinformationmany times,upuntiltheoptimalityconditionsaremet,techniquestoreduce theforwardmodelsimulationcosthavebeenproposed(Jansen,2011; Cardosoetal.,2009;vanDorenetal.,2006).
Multiscale(MS)simulationmethods(Jennyetal.,2003;Houand Wu,1997)havebeenincreasinglyemployedfortheefficientsolution
∗Correspondingauthor.PetrobrasResearchandDevelopmentCenterCENPES,Av.HorácioMacedo950,CidadeUniversitária,RiodeJaneiro,RJ21941-915,
Brazil
E-mail addresses: r.moraes@tudelft.nl (R.J. de Moraes), wesseldezeeuw@hotmail.com (W. de Zeeuw), jrprodrigues@petrobras.com (J.R.P. Rodrigues),
h.hajibeygi@tudelft.nl (H.Hajibeygi),j.d.jansen@tudelft.nl (J.D.Jansen).
ofelliptic(ZhouandTchelepi,2008)andparabolic(Ţeneetal.,2015) equations, more specificallysubsurface flow problemsin highly het-erogeneousporousmedia.Also,manydevelopmentstoextendits ap-plicabilitytoextendedphysicshavebeenobservedintherecentyears (Lieetal.,2017).
Moreover,MSderivativecomputationstrategies,basedonMS for-ward simulation models, has also been subject of study. MS ad-joint formulation (MS-ADJ) for single-phase subsurface flow have been presented in Fu et al. (2010, 2011). MS adjoint computation methods for multiphase flow have also been developed (Krogstad et al., 2011; Moraes et al., 2017). More recently, a mathemati-cal framework for MS computation of derivative information has been developed(deMoraes etal., 2017).Ithas beenhighlightedin
deMoraesetal.(2017);Fuetal.(2010)thatinaccurateMSgradient computationscouldleadtoinaccurategradientdirections.However,as itisindicatedindeMoraesetal.(2017),strategiesthatimprovethe MSforwardsimulationsolution(e.g.refinementoftheMScoarsegrid) resultinbettergradientestimates.Moreover,inFuetal.(2010)itis sug-gestedthataniterativeMSgradientcomputationstrategycouldresolve themultiscalegradientinaccuracies.
https://doi.org/10.1016/j.advwatres.2019.05.016
Received18September2018;Receivedinrevisedform20May2019;Accepted23May2019 Availableonline23May2019
Inthiswork,wedevelopaniterativemultiscalegradient computa-tionstrategywhichconvergestothefine-scalegradientsolution,thus allowingforerrorcontrolandreductionformultiscalegradients.The derivativecomputationmethodisbasedonthegenericmathematical frameworkintroducedindeMoraesetal.(2017).Byaugmentingthe MS modelequation andthe state vectors with the i-MSFV smooth-ingstage,theframeworkiscapable of providingderivative informa-tionatanydesiredaccuracy,uptofine-scaleprecision.The augmen-tationisaddressedbytheimplicitdifferentiationstrategy.Inthe for-mulation,weavoidadditionalcomplexityinvolvedincomputing par-tialderivativesassociatedtothesmoothingstepbyalsoonly approxi-matelysolvingthederivativestateequationassociatedtotheit.Also, the strategy seamlessly accommodates both the Direct and Adjoint methods.
The remaining of this paper is organized as follows. Firstly, we presentanalgebraicformulationforthei-MSFVmethod,suitable for thederivationofthederivativecomputationmethods.Next,wederive theDirectandAdjointmethodstocomputederivativeinformation, fol-lowingthei-MSFVframework,whenthealgorithmsarepresented. Nu-mericalvalidationof themethodagainstnumericaldifferentiationis presented.Thenumericalexperimentsconductedshowthatfine-scale gradientcanbereproducedviathei-MSFVmethodiftheforward simu-lationconvergestoasmallenoughresidualtolerance.Weshow numer-icalevidencethatthereisarelationshipbetweenthegradientquality andthei-MSFVsolutionresidualbycomparingfine-scalegradientand i-MSFVgradientsandrelatingthedifferencebetweenthetwowiththe pressureerrornorm.Concludingremarksarefinallypresentedinthe lastsection.
2. Algebraicandalgorithmicdescriptionofthemultiscale iterativemethod
Weconsiderthesetofequationsthatalgebraicallydescribesthe for-wardsimulationatthefinescale,withoutanyassumptionregardingthe underlyingphysicalmodel,as(deMoraesetal.,2017)
𝐠𝐹(𝐱,θ)=𝟎, (1)
where𝐠𝐹∶ℝ𝑁𝐹×ℝ𝑁𝜃→ ℝ𝑁𝐹 representsthesetofalgebraicforward modelequations,𝐱∈ℝ𝑁𝐹 isthestatevector(which,forsingle-phase flow,containsthegridblockpressures),θ ∈ℝ𝑁𝜃isthevectorof param-eters,andthesubscriptFrefersto‘finescale’.ThereareNFfine-scale cellsandN𝜃parameters.Eq.(1)implicitlyassumesadependencyofthe statevectorxontheparametersθ,i.e.
𝐱=𝐱(θ). (2)
Oncethemodelstateis determined,theobservableresponsesof the forwardmodelarecomputed.Theforwardmodelresponsesmaynot onlydependonthemodelstate,butalsoontheparametersthemselves, andcanbeexpressedas
𝐲𝐹=𝐡𝐹(𝐱,θ), (3)
where 𝐡𝐹∶ℝ𝑁𝐹×ℝ𝑁𝜃→ ℝ𝑁𝑦 represents the output equations (Jansen,2016).ItisassumedthatgFcanbedescribedas
𝐠𝐹(𝐱,θ)=𝐀(θ)𝐱−𝐪(θ), (4) where𝐀(θ)∈ℝ𝑁𝐹×ℝ𝑁𝐹 matrixand𝐪(θ)∈ℝ𝑁𝐹.
Atwo-stagemultiscale(MS)solutionstrategy(Jennyetal., 2003; Wangetal.,2014)canbedevisedbyfirstlycomputingacoarsescale solution
̆𝐠(̆𝐱,θ)=(𝐑𝐀𝐏)̆𝐱−(𝐑𝐪)= ̆𝐀̆𝐱−̆𝐪=̆𝟎, (5)
̆𝐠∶ℝ𝑁𝐶×ℝ𝑁𝜃→ ℝ𝑁𝐶,whereNCisthenumberofcoarsegrid-blocks, andthenanapproximatedfine-scalesolution
𝐠′(𝐱′,̆𝐱,θ)=𝐱′−𝐏̆𝐱=𝟎, (6)
where𝐠′∶ℝ𝑁𝐹×ℝ𝑁𝐶×ℝ𝑁𝜃→ ℝ𝑁𝐹.
Theso-called prolongationoperator𝐏=𝐏(θ)whichis anNF×NC matrixthat maps(interpolates) thecoarse-scalesolutiontothe fine-scale resolution. The so-called restriction operator 𝐑=𝐑(θ) is de-fined asan NC×NF matrixwhich maps thefine scaleto thecoarse scale. Moreinformation about howthese operators are constructed for the Multiscale Finite Volume (MSFV) method can be found in
ZhouandTchelepi(2008);Wangetal.(2014).Let̆𝐱∈ℝ𝑁𝐶bethecoarse scale solution (NC≪NF), and 𝐱′∈ℝ𝑁𝐹 the approximated fine-scale solution.
Theiterativemultiscalestrategy(Hajibeygietal.,2008)canbe de-visedbyconsideringversionsofEqs.(4)and(5)writteninresidualform. Let𝐱𝜈−1beanapproximatesolutiontoEq.(4)atiteration𝜈 −1and
𝐫𝜈−1=𝐪−𝐀𝐱𝜈−1 (7)
bethecorrespondingresidual.Amultiscaleimprovementtothis approx-imationcanbedevisedbywritingEq.(5)inresidualformas
̆𝐠𝜈(̆𝐱𝜈,𝐱𝜈−1,θ)= ̆𝐀̆𝐱𝜈−̆𝐫𝜈−1=̆𝟎, (8) where
̆𝐫𝜈−1=𝐑𝐫𝜈−1, (9)
̆𝐫𝜈−1∈ℝ𝑁𝐶.Here,̆𝐱𝜈isredefinedasthecoarsescalecorrection. Redefin-ingEq.(6),wehave
𝐠′𝜈(𝐱′𝜈,̆𝐱𝜈,θ)=𝐱′𝜈−𝐏̆𝐱𝜈=𝟎, (10)
suchthatx′𝜈 nowrepresentstheapproximatefine-scalecorrectionat iteration𝜈,i.e.,
𝐱𝜈−1∕2=𝐱𝜈−1+𝐱′𝜈 (11)
is anapproximatesolutionofEq.(4)augmentedwiththecorrection fromthecoarse-scalecalculation.
TheapproximatesolutionprovidedbyEq.(11)canbeimprovedif successivesmoothingstepsareemployed(Hajibeygietal.,2008).Let
𝐫𝜈−1
𝜎 =𝐪−𝐀𝐱𝜈−1∕2=𝐪−𝐀 (
𝐱𝜈−1+𝐱′𝜈), (12)
𝐫𝜈−1
𝜎 ∈ℝ𝑁𝐹,bethesmoothedresidualobtainedfromtheapproximation givenbyEq.(11)and
𝐠𝜈
𝜎(𝐱′𝜈,𝐱𝜈𝜎,𝐱𝜈−1,θ)=𝐀𝐱𝜈𝜎−𝐫𝜎𝜈−1=𝟎, (13) a versionof Eq.(4) written in residual form.Here 𝐱𝜈
𝜎∈ℝ𝑁𝐹 is the smoothedfine-scalecorrectionatiteration𝜈.Thesolutionsmoothingis obtainedbysolvingEq.(13)usinganyiterativesolveruptoaprescribed (loose)toleranceor(small)maximumnumberofiterations(Hajibeygi etal.,2008;Ţeneetal.,2015).Thesolutionforagiveniteration𝜈 is,
hence,obtainedfrom
𝐠𝜈 𝑥 ( 𝐱𝜈,𝐱′𝜈,𝐱𝜈 𝜎,𝐱𝜈−1,θ ) =𝐱𝜈−𝐱𝜈−1−𝐱′𝜈−𝐱𝜈𝜎=𝟎, (14) where𝐠𝜈 𝑥∶ℝ𝑁𝐹×ℝ𝑁𝐹×ℝ𝑁𝐹×ℝ𝑁𝐹×ℝ𝑁𝜃→ ℝ𝑁𝐹.
TheMSiterativestrategyisfullydepictedinAlgorithm1,where∥.∥ representsthe2-norm.
InAlgorithm1,𝜖 and𝜖𝜎 are,respectively,theuser-defined toler-ancesfortheouter-loopandsmoothingstep.Thatallowstocontrolthe smoothingstepasarelativeimprovementstartingfromtheMS approx-imatesolution.Aninvestigationofanoptimalrelationshipbetweenthe numberofouterloopsandthenumberofsmoothingstepsispresented inŢeneetal.(2015).
3. Iterativemultiscalegradientcomputation
Forthedevelopmentsthatwillfollowinthissection,itisconvenient towritethesetofequationsthatissolvedineveryiteration,namely
Algorithm1:Iterativemultiscalemethod(Hajibeygietal.,2008).
Input : ̆𝐀,𝐀,𝐑,𝐏,𝐪,𝜖,𝜖𝜎
Output:Approximatesolutionforthelinearsystem𝐀𝐱=𝐪
1 Set𝐱0=𝟎
2 for𝜈 =1,2,…do
3 Compute𝐫𝜈−1=𝐪−𝐀𝐱𝜈−1
4 If∥𝐫𝜈−1∥
∥𝐫0∥ <𝜖,quitwithsolutiongivenby𝐱=𝐱 𝜈−1 5 Solvĕ𝐱𝜈= ̆𝐀−1(𝐑𝐫𝜈−1) 6 Compute𝐱′𝜈=𝐏̆𝐱𝜈 7 Compute𝐫𝜈−1 𝜎 =𝐪−𝐀(𝐱𝜈−1+𝐱′𝜈) 8 Iterativelysolve𝐱𝜈 𝜎=𝐀−1𝐫𝜈−1𝜎 until ∥𝐫𝜎𝜈−1−𝐀𝐱𝜈𝜎∥ ∥𝐫𝜎𝜈−1∥ <𝜖𝜎 9 Update𝐱𝜈=𝐱𝜈−1+𝐱′𝜈+𝐱𝜈 𝜎
Eqs.(8),(10),(13)and(14),inmatrixformas
(15)
Onemustnote,however,that,inthei-MSFVprocedure,Eq.(13)is solvedonlyapproximatelyand,therefore,strictlyspeaking the equa-tionin thethirdrowof Eq.(15) doesnot hold.The ideahere is to describetheprocedureinanalgebraicmanner,ignoringthis approx-imation,inordertofacilitatethepresentationofthederivative calcu-lationalgorithmsinthenextsection.Oncethederivativecalculation methodsareobtainedundertheassumptionthatthealgebraicrelations inEq.(15) hold,thesametypeof smoothingapproach employedin thei-MSFV toresolvehighfrequency errors will be usedin the so-lutionof thederivativeinformation. Thisresults ina practical semi-analyticalalgorithmforderivativecalculationinaniterative formula-tion.Notethatafullyanalyticalprocedurewouldrequirecalculatingthe derivativeofthesmoothingoperatoremployedinthei-MSFV(step8in
Algorithm1),whichcanbequitecomplex,duetoitsnonlinear charac-ter(Frank,2017).Theproposedsemi-analyticalapproachalsobecomes generalandapplicabletoanyiterativeprocedureusedinstep8, regard-lessofitsnature.Atrulyanalyticalderivativemethodwouldrequirethe implementationofderivativecalculationforeachiterativeprocedure used.Moredetailsontheimplicationofthisassumptionarediscussed in3.2.1.
Consideringallequationsthatmustbesolvedforalliterations𝜈 ∈
{1,…,𝑁𝜈},they canbe collated ina super-vector(Rodrigues, 2006; Kraaijevangeretal.,2007;Jansen,2016)fashionas
𝒈(𝒙,θ)= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ̆𝐠1(̆𝐱1,𝐱′0,𝐱0 𝜎,θ ) 𝐠′1(𝐱′1,̆𝐱1,θ) 𝐠1 𝜎 ( 𝐱1 𝜎,𝐱′1,𝐱′0,𝐱0𝜎,θ ) 𝐠1 𝑥 ( 𝐱1,𝐱1 𝜎,𝐱′1,𝐱0,θ ) ̆𝐠2(̆𝐱2,𝐱′1,𝐱1 𝜎,θ ) 𝐠′2(𝐱′2,̆𝐱2,θ) 𝐠2 𝜎 ( 𝐱2 𝜎,𝐱′2,𝐱′1,𝐱1𝜎,θ ) 𝐠2 𝑥 ( 𝐱2,𝐱2 𝜎,𝐱′2,𝐱1,θ ) ⋮ ̆𝐠𝑁𝜈(̆𝐱𝑁𝜈,𝐱′𝑁𝜈−1,𝐱𝑁𝜈−1 𝜎 ,θ ) 𝐠′𝑁𝜈(𝐱′𝑁𝜈,̆𝐱𝑁𝜈,θ ) 𝐠𝑁𝜈 𝜎 ( 𝐱𝑁𝜈 𝜎 ,𝐱′𝑁𝜈,𝐱′𝑁𝜈−1,𝐱𝑁𝜎𝜈−1,θ ) 𝐠𝜈 𝑥 ( 𝐱𝑁𝜈,𝐱𝑁𝜈 𝜎 ,𝐱′𝑁𝜈,𝐱𝑁𝜈−1,θ ) ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ =𝟎, (16)
andthesuper-statevectordefinedas
𝒙(θ)= ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ ̆𝐱1 𝐱′1 𝐱1 𝜎 𝐱1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ 𝑇 … ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ ̆𝐱𝑁𝜈−1 𝐱′𝑁𝜈−1 𝐱𝑁𝜈−1 𝜎 𝐱𝑁𝜈−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ 𝑇 ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ ̆𝐱𝑁𝜈 𝐱′𝑁𝜈 𝐱𝑁𝜈 𝜎 𝐱𝑁𝜈 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ 𝑇⎞ ⎟ ⎟ ⎟ ⎟ ⎠ 𝑇 . (17)
Notethatweusebold-italicinthenotationtorepresentsuper-vectors. It is discussed in Rodrigues (2006); de Moraes et al. (2017),
deMoraesetal.howanyderivativeinformationcanbeefficiently com-putedfromthepreandpostmultiplicationofthesensitivitymatrixG,
𝐆∈ℝ𝑁𝑦×𝑁𝜃,byarbitrarymatricesas 𝐖𝐆𝐕=−𝐖𝜕𝐡 𝜕𝒙 ( 𝜕𝒈 𝜕𝒙 )−1 𝜕𝒈 𝜕θ𝐕+𝐖𝜕𝐡𝜕θ𝐕, (18)
where V (of order 𝑁θ×𝑝) and W (of order m×Ny) are arbitrary matrices, defined based on thederivative information one wants to obtain.
Eq. (18) requires the partial derivative of the model equations withrespecttotheparameters.FromEqs.(5)and(6),asdiscussedin
deMoraesetal.(2017),itfollowsthat
𝜕̆𝐠𝜈 𝜕θ =[ 𝜕𝐑𝜕θ(𝐀𝐏)+𝐑𝜕𝐀𝜕θ𝐏+(𝐑𝐀)𝜕𝐏𝜕θ ] ̆𝐱𝜈−𝜕𝐑 𝜕θ𝐪+ −𝐑𝜕𝐪 𝜕θ+( 𝜕𝐑𝜕θ𝐀+𝐑𝜕𝐀𝜕θ ) 𝐱𝜈−1, (19) 𝜕𝐠′ 𝜕θ 𝜈 =−𝜕𝐏 𝜕θ̆𝐱𝜈, (20)
andderivingEq.(13)withrespecttoθ
𝜕𝐠𝜈 𝜎
𝜕θ = 𝜕𝐀𝜕θ𝐱𝜎𝜈+𝜕𝐀𝜕θ𝐱′𝜈−𝜕𝐪𝜕θ+𝜕𝐀𝜕θ𝐱𝜈−1. (21) Also,fromEq.(14),itfollowsthat
𝜕𝐠𝜈 𝑥
𝜕θ =𝟎 (22)
ForthesakeofsimplicityandincoherencewiththeMSFVmethod used inthenumericalexperiments,thedependencyof Ron θis ne-glected.
Now,theordertheoperationsinEq.(18)areevaluateddefinethe Direct andAdjoint methods.Thederivation ofboth methodswillbe discussedinthenextsections.
3.1. Thedirectmethod
IfWisfactoredoutinEq.(18),itcanberewrittenas
𝐆𝐕= 𝜕𝐡 𝜕𝒙𝐙+𝜕𝐡𝜕θ𝐕, (23) where 𝐙=− ( 𝜕𝒈 𝜕𝒙 )−1 𝜕𝒈 𝜕θ𝐕, (24) issolvedfrom ( 𝜕𝒈 𝜕𝒙 ) 𝐙=−𝜕𝒈 𝜕θ𝐕. (25)
ThelinearsystemdescribedinEq.(25)canbere-writtenina block-wisefashionforeachiteration𝜈:
(26) where,fromEqs.(8),(10),(13),and(14)
𝜕𝐠𝜈 𝜕𝒙𝜈 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝜕̆𝐠𝜈 𝜕𝒙𝜈 𝜕𝐠′𝜈 𝜕𝒙𝜈 𝜕𝐠𝜈 𝜎 𝜕𝒙𝜈 𝜕𝐠𝜈 𝑥 𝜕𝒙𝜈 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ ̆𝐀 𝟎 𝟎 𝟎 −𝐏 𝐈 𝟎 𝟎 𝟎 𝐀 𝐀 𝟎 𝟎 −𝐈 −𝐈 𝐈 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ , (27) and 𝜕𝐠𝜈 𝜕𝒙𝜈−1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 𝜕̆𝐠𝜈 𝜕𝒙𝜈−1 𝜕𝐠′𝜈 𝜕𝒙𝜈−1 𝜕𝐠𝜈 𝜎 𝜕𝒙𝜈−1 𝜕𝐠𝜈 𝑥 𝜕𝒙𝜈−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 𝟎 𝟎 𝟎 𝐑𝐀 𝟎 𝟎 𝟎 𝟎 𝟎 𝟎 𝟎 𝐀 𝟎 𝟎 𝟎 −𝐈 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ . (28)
Thepartitioninglinesindicatewhichmatrixandvectortermsbelongto eachiteration.
SubstitutingEq.(27)and(28)inEq.(26),followsthat
(29) Foreveryiteration,thelinearsystemthatmustbesolvedforthe coarse-scaleequationis ̆𝐙𝜈= ̆𝐀−1 ( −𝜕̆𝐠 𝜈 𝜕θ𝐕−𝐑𝐀𝐙𝜈−1𝑥 ) , (30)
forthefine-scaleapproximatesolutionequation
𝐙′𝜈=−𝐏̆𝐙𝜈−𝜕𝐠′𝜈
andforthesmoothingequation 𝐙𝜈 𝜎=𝐀−1 ( −𝜕𝐠 𝜈 𝜎 𝜕θ𝐕−𝐀𝐙′𝜈−𝐀𝐙𝜈−1𝑥 ) . (32)
Aspointedoutabove,itwouldnotbefeasibletofullysolveEq.(32). In order to obtain a practical derivative calculation method, the same smoothingprocedure used in the i-MSFVis applied here,i.e.,
𝐙𝜈
𝜎 is obtained by solving Eq. (32) using any iterative solver up to a prescribed (loose) tolerance or a (small) maximum number of iterations.
Finally,
𝐙𝜈
𝑥=𝐙𝜈−1𝑥 +𝐙′𝜈+𝐙𝜈𝜎. (33)
Observingthathonlydependson𝐱=𝐱𝑁𝜈,Eq.(23)canbesimplified to
𝐆𝐕= 𝜕𝐡
𝜕𝐱𝐙𝑁𝑥𝜈+𝜕𝐡𝜕θ𝐕. (34) NowtheDirectmethodfortheiterativemultiscalegradient compu-tationcanbefullydetermined.ItisdepictedinAlgorithm2.Notethat referencestosuperscript−1correspondtozeroterms.
Algorithm2:PostmultiplicationofGbyVviathedirectmethod.
Input :𝐑,𝐀,𝐏,𝒙, 𝜕𝐀
𝜕θ, 𝜕𝐡𝜕𝐱,𝜕𝐡𝜕θ,𝐕,𝜖𝜎
Output:The𝐆𝐕product
1 foreach𝑗=1,2,…,𝑛do 2 foreach 𝜈 =0,1,…,𝑁𝜈do 3 Computeβ =𝜕𝐏 𝜕θ̆𝐱𝜈𝐕.,𝑗;// Algorithm 3, in [19] where 𝐦=𝐕.,𝑗 4 Computeα =𝐑𝐀β +𝐑 ( 𝜕𝐀 𝜕θ𝐱𝜈−1−𝜕𝐪𝜕θ+𝜕𝐀𝜕θ𝐏̆𝐱𝜈 ) 𝐕.,𝑗 5 Solve ̆𝐙𝜈.,𝑗= ̆𝐀−1(−α −𝐑𝐀𝐙 𝑥𝜈−1.,𝑗 ) 6 Compute𝐙′𝜈 .,𝑗=−𝐏̆𝐙𝜈.,𝑗+β 7 Compute𝜕𝐠 𝜈 𝜎 𝜕θ =𝜕𝐀𝜕θ𝐱𝜈𝜎+𝜕𝐀𝜕θ𝐱′𝜈−𝜕𝐪𝜕θ+𝜕𝐀𝜕θ𝐱𝜈−1 8 Computeδ = ( −𝜕𝐠 𝜈 𝜎 𝜕θ𝐕.,𝑗−𝐀𝐙′𝜈.,𝑗−𝐀𝐙𝑥𝜈−1.,𝑗 )
9 Iterativelysolve𝐙𝜎𝜈.,𝑗=𝐀−1δuntil∥δ −𝐀𝐙𝜎 𝜈 .,𝑗∥ ∥δ ∥ <𝜖𝜎 10 Compute𝐙𝑥𝜈.,𝑗=𝐙𝑥.,𝑗𝜈−1+𝐙′𝜈.,𝑗+𝐙𝜎𝜈.,𝑗 11 Compute(𝐆𝐕).,𝑗= 𝜕𝐡 𝜕𝐱𝐙𝑥𝑁.,𝑗𝜈+𝜕𝐡𝜕θ𝐕.,𝑗
3.2. Theadjointmethod
Now,ifVisfactoredoutinEq.(18),itcanberewrittenas
𝐖𝐆=𝐙𝜕𝒈 𝜕θ+𝐖𝜕𝐡𝜕θ, (35) where 𝐙=−𝐖𝜕𝐡 𝜕𝒙 ( 𝜕𝒈 𝜕𝒙 )−1 (36) issolvedfrom 𝐙 ( 𝜕𝒈 𝜕𝒙 ) =−𝐖𝜕𝐡 𝜕𝒙. (37)
ThelinearsystemdescribedinEq.(37)canbere-writtenina block-wisefashionforeachiteration𝜈 as
(38)
ThestructureofEq.(38)showsthatitcanbesolvedviaback substi-tution,i.e.thesolutionofthegradientinformationisbackwardinthe iterations.
SubstitutingEqs.(27)and(28)inEq.(38),followsthat
(39) From Eq. (39) and recalling that h only depends on 𝐱=𝐱𝑁𝜈, the (transposed) linear system that must be solved to compute Z𝜈 is
givenby
(40)
for𝜈 ≠ N𝜈,while𝐙𝑁𝜈iscalculatedfrom ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ ̆𝐀𝑇 −𝐏𝑇 𝟎 𝟎 𝟎 𝐈 𝐀𝑇 −𝐈 𝟎 𝟎 𝐀𝑇 −𝐈 𝟎 𝟎 𝟎 𝐈 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ( ̆𝐙𝑁𝜈)𝑇 ( 𝐙′𝑁𝜈)𝑇 ( 𝐙𝑁𝜈 𝜎 )𝑇 ( 𝐙𝑁𝜈 𝑥 )𝑇 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 𝟎 𝟎 𝟎 −( 𝜕𝐡 𝜕𝐱 )𝑇 𝐖𝑇 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ . (41)
Theequationthatcomputesthe”adjointstate” associatedwith𝐠𝜈𝑥reads ( 𝐙𝜈 𝑥 )𝑇=(𝐙𝜈−1 𝑥 )𝑇 −𝐀𝑇𝐑𝑇(̆𝐙𝜈−1)𝑇−𝐀𝑇(𝐙𝜈−1𝜎 )𝑇,𝜈 ≠ 𝑁𝜈, (42) and ( 𝐙𝑁𝜈 𝑥 )𝑇 =−( 𝜕𝐡 𝜕𝐱 )𝑇 𝐖𝑇. (43)
The“adjointstate” for𝐠𝜈
𝜎iscalculatedfrom (
𝐙𝜈
𝜎)𝑇=𝐀−𝑇(𝐙𝜈𝑥)𝑇. (44)
Asdiscussedforthedirectmethod,intheproposedderivative cal-culationmethod,Eq.(44)issolvedasasmoothingsteponly,usingany iterativesolveruptoaprescribedtoleranceornumberofiterations.
Forg′𝜈,theadjointstateisgivenby (
𝐙′𝜈)𝑇=(𝐙𝜈
𝑥)𝑇−𝐀𝑇(𝐙𝜈𝜎)𝑇, (45)
andfinallyfor̆𝐠𝜈
̆𝐙𝜈=(̆𝐀−𝑇)(𝐏𝑇𝐙′𝜈). (46) Eq.(35)canbewrittenasasumwhereeachtermcorrespondstothe contributionofoneiteration
𝐖𝐆= 𝑁𝜈 ∑ 𝜈=0 ( ̆𝐙𝜈𝜕̆𝐠𝜈 𝜕θ +𝐙′𝜈 𝜕𝐠′𝜈 𝜕θ +𝐙𝜈𝜎 𝜕𝐠𝜈 𝜎 𝜕θ ) +𝐖𝜕𝐡 𝜕θ. (47)
The algorithm can now be fully defined and is presented in 3. Note the use of the notation (WG)𝜈 to denote the partial sums in
Eq.(47) andthatreferencestosuperscript𝑁𝜈+1correspondtozero terms.
Algorithm3:PremultiplicationofGbyWviatheadjointmethod.
Input :𝐑,𝐀,𝐏,𝒙,𝜕𝐀
𝜕θ,𝜕𝐡𝜕𝐱, 𝜕𝐡𝜕θ,𝐖,𝜖𝜎
Output:The𝐖𝐆product
1 foreach𝑖=1,2,…,𝑚do 2 foreach 𝜈 =𝑁𝜈,…,1,0do 3 Set(𝐙𝜈 𝑥 )𝑇 𝑖,.= ⎧ ⎪ ⎨ ⎪ ⎩ −( 𝜕𝐡 𝜕𝐱 )𝑇 𝐖𝑇,if𝜈 =𝑁 𝜈 ( 𝐙𝜈+1 𝑥 )𝑇 𝑖,.−𝐀𝑇 ( 𝐙𝜈+1 𝜎 )𝑇 𝑖,.−𝐀𝑇𝐑𝑇 ( ̆𝐙𝜈+1)𝑇 𝑖,.,otherwise 4 Iterativelysolve(𝐙𝜈 𝜎 )𝑇 𝑖,.=𝐀−𝑇 ( 𝐙𝜈 𝑥 )𝑇 𝑖,.until ∥𝐙𝜈𝑥−𝐀𝐙𝜎𝜈.,𝑗∥ ∥𝐙𝜈𝑥∥ <𝜖𝜎 5 Compute(𝐙′𝜈)𝑇 𝑖,.= ( 𝐙𝜈 𝑥 )𝑇 𝑖,.−𝐀𝑇 ( 𝐙𝜈 𝜎 )𝑇 𝑖,. 6 Solve ̆𝐙𝜈𝑖,.=(̆𝐀−𝑇)(𝐏̆𝐙𝜈+1 𝑖,. ) 7 Computeα𝑇= ̆𝐙𝜈𝑖,.(𝐑𝐀)−𝐙′𝜈 𝑖,. 8 Computeβ =α 𝜕𝐏 𝜕θ̆𝐱𝜈;// Algorithm 4 in [19] where 𝐦𝑇 =α𝑇 9 Computeγ =𝐑𝜕𝐀 𝜕θ𝐏̆𝐱𝜈−𝐑𝜕𝐪𝜕θ+𝐑𝜕𝐀𝜕θ𝐱𝜈−1 10 Compute𝜕𝐠 𝜈 𝜎 𝜕θ = 𝜕𝐀𝜕θ𝐱𝜎𝜈+𝜕𝐀𝜕θ𝐱′𝜈−𝜕𝐪𝜕θ+𝜕𝐀𝜕θ𝐱𝜈−1 11 Update(𝐖𝐆)𝜈𝑖,.=(𝐖𝐆)𝜈+1𝑖,. +β +̆𝐙𝜈𝑖,.γ +𝐙𝜈𝜎𝑖,.𝜕𝐠 𝜈 𝜎 𝜕θ 12 Update𝐖𝐆=(𝐖𝐆)0+𝐖𝜕𝐡 𝜕θ
3.2.1. Remarksabouttheframework
Analternativeformulationforthei-MSFVformulationhasbeen pre-viouslyproposedandinvestigatedinFrank(2017).Inthatwork,both the state and themodel equation vectors explicitly account for the smoothingstage.Theformulationhereproposedis basedon two ob-servations.Firstly,theimplementationoftheaforementionedvariant, although offersmore controloverthegradient quality,relies on the abilityofcomputingpartialderivativematricesofsmoothingstepwith respecttotheparameters.Morespecifically,itimpliesontheknowledge ofhowthepreconditionMofthesystemmatrixAisbuilt.Forsimpler iterativestrategies,e.g.Jacobi,whichconstructionofMcanbesimple, thecomputationof 𝜕𝐌
𝜕θ isrelativelysimple.However,simpleriterative
methodsareusuallylessefficient.Also,therequirementofknowingthe constructionof Mhamperstheutilization of‘black-box’ typeof pre-conditioners.Secondly,ithasbeenshown inŢene etal.(2015)that onlyalimitednumberofsmoothingstepsarenecessarytoresultinan efficienti-MSFVsolutionstrategy.Hence,notmuchextracontrolwould beachieved.
Notethatthelinearsolvers employedinlines9and4of, respec-tively,Algorithms2and3,arethesamesolversemployedinthe solu-tionoftheforwardsimulation,usingthesameprescribed(loose) toler-ance.Hence,thealgorithmssharethesamecomputationaladvantages bysolvingfortheapproximatedderivativeinformationarisingfromthe smoothingstep.
Thebackwardalgorithmrequiresstoringallintermediatestates gen-eratedduringtheiterationsintheforwardrun.Italsorequiressolving manysystemsofequationsforeachbackwardtime-step.Iftheiteration processintheforwardrungoesuntilmachineprecisionisreached,then essentiallythefullycoupledsystemhasbeensolvedtofine-scale preci-sionanditmightbemorebeneficialtoneglecttheiterationhistoryand aimtosolvethefine-scalesystemofadjointequationsinamoreefficient way,giventhatthederivativecomputationproblemislinear.However, wehighlightthat,asithasbeenshownintheliterature(Fonsecaetal., 2015;deMoraesetal.,2017),approximategradientscomputedfrom
approximatesolutionsarealreadysufficienttoefficiently/successfully leadtheoptimizationprocesstotheoptimalsolution.Howaccuratethis gradientwillneedtobeisdependentondifferentaspects,e.g.the opti-mizationalgorithm,howearly/latetheoptimizationprocessisinterms offindingthemaximum/minimumofagivenobjectivefunction,among others.Thealgorithmhereproposedhastheadvantageofcontrolling thegradientquality,aswillbeshowninthenumericalexperiments pre-sentedinthenextsection.
4. Numericalexperiments
Giventhefundamentalnatureofthedevelopmentspresentedhere, wefocusthenumericalexperimentsonthevalidationofthe computa-tionandassessmentofthegradientqualityprovidedbytheiterative multiscalemethod.
Ourexperimentswillbebasedontheevaluationofthegradientof amisfitobjectivefunctionwithnoregularizationterm(Oliveretal., 2008) 𝑂(θ)=1 2 ( 𝐡(𝐱,θ)−𝐝𝑜𝑏𝑠)𝑇𝐂𝐷−1(𝐡(𝐱,θ)−𝐝𝑜𝑏𝑠), (48) withagradient ∇𝛉𝑂=𝐆𝑇𝐂−1 𝐷(𝐡(𝐱,θ)−𝐝𝑜𝑏𝑠), (49)
where𝐂𝐷∈ℝ𝑁𝑌×𝑁𝑌 istheso-calleddatacovariancematrix.Here,CD willbeconsideredadiagonalmatrixgivenbyOliveretal.(2008)
𝐂𝐷=𝜎2𝐈, (50)
where𝜎2isthevarianceofthedatameasurementerror.
Inallexperiments,thefittingparametersarecell-centered perme-abilities.Theobservedquantity,dobs,isthefinescalepressureatthe locationof(non-flowing)observationwells,therefore
𝜕𝐡
𝜕𝐱′ =𝐈, (51)
and
𝜕𝐡
𝜕̆𝐱 =𝟎. (52)
In all experiments,the standard deviation of the pressure measure-menterroris 𝜎 ≈ 0.03(notethatthemeasurementerror isalso non-dimensional).Thisrepresentsa(veryaccurate)measurementerrorin therangeofthoseusuallyemployedinsyntheticstudycases(seee.g.
Oliveretal.,2008).
Notethatthewellsarecontrolledbybottom-holepressure,expressed intermsofanon-dimensionalpressure,i.e.,
𝑝𝐷= 𝑝 −𝑝𝑝𝑟𝑜𝑑
𝑝𝑖𝑛𝑗−𝑝𝑝𝑟𝑜𝑑, (53)
wherepinjandpprodaretheinjectionandproductionpressures, respec-tively.Inalltheexperiments,𝑝𝑖𝑛𝑗=1.0and𝑝𝑝𝑟𝑜𝑑=0.0,thegrid-block dimensionsareΔ𝑥=Δ𝑦=Δ𝑧=1mandthefluidviscosityis1.0× 10−3 Pas.Inaddition,inallthefollowingtestcases,wellbasisfunctionsare included.
Themetricutilizedtoassessthei-MSFVgradientqualityistheangle betweenfine-scaleandi-MSFVnormalizedgradients,i.e.,
𝛼 =cos−1(∇𝑇θ̂𝑂𝐹 𝑆∇θ̂𝑂𝑀𝑆). (54) Here, ∇θ̂𝑂𝐹 𝑆= ∇θ𝑂𝐹 𝑆 ‖‖∇θ𝑂𝐹 𝑆‖‖2 (55) and ∇θ̂𝑂𝑀𝑆= ∇θ𝑂𝑀𝑆 ‖‖∇θ𝑂𝑀𝑆‖‖2 . (56)
Also,∇θ𝑂𝐹 𝑆 and∇θ𝑂𝑀𝑆 denotethefine-scaleandMSanalytical gradients,respectively.Asaminimumrequirement,acceptableMS gra-dientsareobtainedif𝛼 ismuchsmallerthan90o(Fonsecaetal.,2015).
Andtoproveourhypothesis,weparticularlyinterestedin observing thebehaviourofthemetricasmoreaccuratei-MSFVsolutionsare com-puted.
Inouri-MSFVimplementation,theiterativeprocessiscontrolledby theouterloopresidual𝜖 andthepre-conditionersmoothererror toler-ance 𝜖𝜎.TheKrylovsubspacebiconjugategradientstabilizedmethod (BiCGSTAB,Saad,2003)isemployedinthesmoothingstage.
4.1. Validationexperiments
Beforefocusinginthevalidation,inthissectionwewillvalidatethe iMS-gradientmethodagainstthenumericaldifferentiationmethodwith ahigherorder,two-sidedTaylorapproximation
∇θ𝑂𝑖= 2𝛿𝜃1 𝑖 ( 𝑂(𝜃1,⋯,𝜃𝑖−1,𝜃𝑖+𝛿𝜃𝑖,𝜃𝑖+1,⋯,𝜃𝑁𝜃) −𝑂(𝜃1,⋯,𝜃𝑖−1,𝜃𝑖−𝛿𝜃𝑖,𝜃𝑖+1,⋯,𝜃𝑁𝜃) ) (57) whereweconsider𝛿 tobeamultiplicativeparameterperturbation.We definetherelativeerroras
𝜀=||∇θ𝑂𝑁𝑈𝑀−∇θ𝑂𝑖𝑀𝑆||2 ||∇θ𝑂𝑖𝑀𝑆||2
(58) Here ∇θ𝑂𝑁𝑈𝑀 is obtainedbyperformingtheproperamountof iter-ative multiscale reservoirsimulations requiredto evaluateEq.[57]. ∇θ𝑂𝑖𝑀𝑆 isobtainedbyemployingtheiterativeDirectorAdjoint
gra-dientmethod.
Toevaluatethecorrectnessoftheproposediterativegradient com-putationmethods,aswellastheirimplementation,weinvestigatethe lineardecreaseoftherelativeerror𝜀bydecreasingtheparameter per-turbation𝛿 from10−1to10−4.Thisinvestigationiscarriedoutintwo
distinctexamples.Bothhaveafinegridof21×21gridblocks.We em-ploya7×7coarseningratio,givinga3×3coarsegrid.Thereference twin-experimentisgeneratedwithpermeabilityrealizationnumber992.
Fig.1illustratesthefine-,coarse-anddual-gridcellsalongwiththe ref-erencepermeability.
Nextwedeterminethewellpositions.Weusetheso-calledquarter wellspot.Here,twoobservationwellsareplacednearoperatingwells. ThefullspecificationscanbefoundinTable1.
The results of this experimentarefoundin Fig.2. Here,we use a singleouter iteration.Weuse a verytight smoothingtoleranceof
𝜖𝜎=5× 10−8 toensurethatthenumericalgradient methodproduces
accurateenoughgradients.Firstofallwecanseethatthefine-scale nu-mericalgradientmethodandtheiterativeMS-gradientmethodsareof thesameorderofaccuracywithrespecttotheperturbation𝛿,forall differentcasesconsidered.Inallexperiments,wecanseethelinear de-creasingbehaviouroftherelativeerrorvalues𝜀astheperturbation𝛿
decreases.Also,notethattheAdjointandDirectmethodsprovidethe analyticalgradientwiththesamelevelofaccuracy.Inthefirst exper-iment,ahomogeneoustestcasewiththepermeabilityvalueof1.0𝑒−13
isused.Thesecondexperimentindicatesthecorrectnessofthemethod whenitisappliedtoheterogeneousporousmediaproblems.
4.2. Investigationofi-MSFVconvergencebehaviourandgradientquality
Inthis investigation,theerror metricgivenbyEq. (54)is evalu-atedforawhole ensembleof heterogeneouspermeability fields.The
Table1
Validationexperimentswellsetup.
Well Fine scale position (I, J) Well Type INJE (1, 1) Injection PROD (21, 21) Production OBS1 (3, 3) Observation OBS2 (19, 19) Observation
Fig.1. Validationexperimentssetup.(a)Primal(boldblack line),dual(identifiedbythebluecells)andfinegrids(thin blacklines).(b)Referencepermeabilityfieldusedinthetwin experiment.(Forinterpretationofthereferencestocolourin thisfigurelegend,thereaderisreferredtothewebversionof thisarticle.)
Fig.2. Validationofthei-MSFVgradientcomputationmethodcomparedwithanumericalgradient.(a)representsthehomogeneoustestcase,while(b)represents theheterogeneoustestcase.
Fig.3. Fourdifferentpermeabilityrealizationsfromtheensembleof1000membersusedinthe2Dnumericalexperiments.
ensembleisgeneratedviathedecompositionofareference permeabil-ityimageusingPrincipalComponentAnalysisparameterization.Fig.3
illustrates4differentpermeabilityrealizationsfromtheensemble.See
Jansen(2013)formoredetails.
Thefine-scaleandcoarsegridscontain21x21and7x7cells, re-spectively.Aninvertedfive-spotwellpatternisemployed,whilefour observationwellsareplaced closetotheproduction wells.The well configurationisdepictedinTable2.
Thesyntheticobservedpressuresattheobservationwelllocations arecreatedviatheclassicaltwin-experimentstrategy,wherethe per-meabilityfieldisextractedfromthe1000geologicalrealizations(the firstoneinFig.3).
TherobustnessofthemethodisillustratedinFig.4.Itispossible toobservethattheanglebetween fine-scalegradientandthei-MSFV gradientissmallerthetighterwemaketheouter-loopresidualtolerance. Moreover,thevariancealsogoestoalmostzeroaswesettheresidual toleranceto1𝑒−4.Forthissetofrelativelysmallpermeabilitycontrast,
perfectalignmentwithfine-scalegradientisreachedifthetoleranceis
Table2
Well configuration for the homogeneous, two-dimensionalcase.
Well Fine scale position (I, J) Well type INJE (11, 11) Injection PROD1 (1, 1) Production PROD2 (21, 1) Production PROD3 (1, 21) Production PROD4 (21, 21) Production OBSWELL1 (3, 3) Observation OBSWELL2 (19, 3) Observation OBSWELL3 (3, 19) Observation OBSWELL4 (19, 19) Observation
setto1𝑒−5.Wehighlightthatthepermeabilitycontrastofthisensemble
isnothigh.Next,weassesstherobustnessforthemethodforgeological settingswithhigherpermeabilitycontrasts.
MSFV 1e-1 1e-2 1e-3 1e-4 1e-5 1e-6
0
10
20
30
(-)
α
(
o)
Fig.4.Box-plotillustratingtheangle𝛼 betweenfine-scalegradientandi-MSFV gradientcomputedforthe1000memberensembleasafunctionofthe outer-looptoleranceerror𝜖.“Noiteration” isequivalenttotheMSFVgradient com-putationpresentedinde Moraes et al. (2017) .
4.3. Robustnesswithrespecttoheterogeneitycontrastanddistribution
Inorder tofurtherexplore thepointabouttherobustnessof the methodwithrespecttoheterogeneitycontrastanddistribution,foursets of20equiprobablerealizationsoflog-normallydistributedpermeability fieldswithasphericalvariogramanddimensionlesscorrelationlengths ofΨ1=0.5andΨ2=0.02aregeneratedusingsequentialGaussian
sim-ulations(Remyetal.,2009).Foreachset,thevarianceandthemean ofln(k)are2.0and3.0,respectively,wherekisthegridblock perme-ability.AsdepictedinFig.5,fortherealizationswithalongcorrelation length,theanglesbetweenthepermeabilitylayersandthehorizontal axisare0o,15o,and45o.Apatchy(smallcorrelationlength)patternis alsoconsidered(Fig.5d).Comparedwiththepreviousset,the perme-abilitycontrastismuchhigherinthiscase.
Thefine-scaleandcoarsegridscontain100×100and20×20cells, respectively.Thewellconfigurationutilizedinthisnumerical experi-mentisdepictedinTable3.
Theobserveddataisgeneratedfromatwin-experimentassociated with(thefirst)permeabilityrealizationofeachset.
Inthisexperiment,𝜖 =1.0𝑒−6and𝜖
𝜎=1.0𝑒−1.Thebox-plotshown inFig.6summarizestherequiredtotalnumberofsmoothingsteps,for allouteri-MSFVsteps.
Thegridorientationeffect(Azizetal.,1993)impactonthe perfor-manceof thei-MSFVmethodisclearin thisexample.Themorethe heterogeneityorientationisalignedwiththefloworientation,theless
Table3
WellconfigurationforCase2.
Well Fine scale position (I, J) Well type INJE (1, 1) Injection PROD (100, 100) Production OBSWELL1 (3, 3) Observation OBSWELL2 (98, 98) Observation
Fig.6. Box-plotrepresentingthetotalnumberofsmoothingiterationsrequired tocomputethemisfitOFgradientforthedifferentpermeabilityensembleswith correlationangles0o,15oand45oandwithasmallcorrelationlength(patchy).
isthenumberofrequirediterations.Also,inrelationtoFig.7,themore challenging theforwardproblem,themore challengingitis to com-putei-MSFVgradientinaccordancetothefine-scalegradient. Never-theless,inallcases,almostalli-MSFVrealizationgradientsareperfectly alignedwithfine-scalegradient, demonstratingtherobustnessof the method.
4.4. SPE-10comparativetestcase
Now, weinvestigate the performanceof ourmethod in the SPE-10 comparative case (Christie et al., 2001), regarded as challeng-ingmodelfor upscaling(Durlofsky,2005)andmultiscalesimulation (Hajibeygi,2011)techniques.Here,weconsiderthe2-Dflow simula-tionofbothtopandbottomlayeroftheoriginal3Dmodel,which per-meabilityfieldsareillustratedinFig.8.
Thefinegriddimensionsis60×220,whileweemploya12x20 coarse grid in theMS simulation.A quarterfive-spotwell setting is
Fig.5. Permeabilitydistributionoffourdifferentrealizationstakenfromthesetsof20geostatisticallyequiprobablepermeabilityfieldswithdifferentcorrelation angles(a-c).Also,apatchyfield(d)withasmallcorrelationlengthisconsidered.
Fig.7. Illustrationofgradientqualityimprovementwhenan i-MSFVgradientcomputationstrategyisemployedin compar-isontoaMSFVcomputationstrategy.Thex-axisrepresentthe angle𝛼 betweenfine-scaleandtheMSFVgradient(illustrated bybluecrosses)andthei-MSFVgradientwitherrortolerance
𝜖 =10−6(illustratedbyorangecircles).(Forinterpretationof
thereferencestocolourinthisfigurelegend,thereaderis re-ferredtothewebversionofthisarticle.)
Fig.8. SPE-10comparativetestcase:top(a)and bot-tom(b)layerpermeabilityfields.
considered.Fourobservationwellsaredeliberatelypositionedinlow permeabilityregions,surroundedbyhighpermeabilityregions.Thewell positionsaredescribedinTable4.
Weemphasizethatithasbeenreported(seee.g.Hajibeygi,2011) thattheMSFVmethodprovidesnon-monotonepressuresolutionswhen solvingthepressureforthebottomlayer(Fig.8b).Strategiestoimprove theMSFVsolution,amongthemthei-MSFV(HajibeygiandJenny,2011) hereconsidered,arenecessarytoaddresstheissue.
Inthistwinexperiment,thereferencepermeabilityfieldusedto com-putetheobservationsisconsideredhomogeneouswithavalueof1e-13.
Table4
WellconfigurationfortheSPE-10comparativetest case.
Well Fine scale position (I, J) Well type INJE (1, 220) Injection PROD (60, 1) Production OBSWELL1 (33, 5) Observation OBSWELL2 (28, 50) Observation OBSWELL3 (28, 83) Observation OBSWELL4 (43, 204) Observation
Fig.9. i-MSFVgradientquality(angle𝛼 betweenfine-scaleandi-MSFVgradients)asafunctionofresidualerror𝜖 fortheSPE-10toplayer(a)andbottomlayer(b).
In this experiment, we fix 𝜖𝜎=1.0𝑒−2 and vary 𝜖 =1.0𝑒−2, 1.0𝑒−4,1.0𝑒−6 to evaluate whether the same convergence behaviour
of the i-MSFV gradient toward fine-scale gradient direction is observed.
Onceagainweobservethatthemethodprovidesaccurategradients, evenconsideringthischallenginggeologicalsetting,forboththetop andbottomlayersofthemodel(seeFig.9).
4.5. Discussion
Basedon theresults acquired from thedifferent numerical cases ofincreasingcomplexitywedemonstratedthatournewlyintroduced method can provideaccurate gradients,up tofine-scale accuracy,if smallenoughresidualtolerancesareemployedinthei-MSFVforward simulation.Also,itisshown thatonlyafewi-MSFV/smoothing iter-ationsarenecessarytohavereasonablyaccurategradients.However, wehighlightthat,fromanoptimizationpointofview,gradientsthat arefullyinaccordance withfine-scalegradientarenotnecessary.As longas thegradient roughlypoints tothecorrect up/downward di-rection,theoptimizationprocessisabletoprogresstowardsthe maxi-mum/minimum.Thisisparticularlytrueintheearlyiterations,where thereisamoresteepregionoftheobjectivefunction.Ontheotherhand, astheoptimizationprocessapproachestheoptimum,moreaccurate gra-dientsarerequiredforaprecisestoppingcriteriaevaluation.The con-vergentbehaviourofthei-MSFVgradient,whichisdemonstratedtobe directlyrelatedtotheouter-residualtolerance,whichisanestimate de-fineda-priori,shouldallowforanerrorcontrolofthegradient computa-tion,whichultimatelywouldallowcontroloftheoptimizationprocess performance.
Thegradientcomputationmethodshere,bydesign,directlyinherits thecharacteristicsofthei-MSFVintroducedinHajibeygietal.(2008). Thesolution of thebackward(transposed)system of equations mir-rorsthesolutionstrategyofthelinearsystemofequationsthatarises from the forward problem. Even the system matrices from the for-wardandthebackwardsystem ofequations share thesame proper-ties(e.g.conditioningnumber).Therefore,in principle,allstudiesin theliterature relatedtothe robustnessandefficiencyof thei-MSFV methodwithrespecttogridsize,levelofheterogeneity,amongothers, shouldalsobevalidinthecontextofthegradientcomputationmethods hereproposed.Nevertheless, furthernumericalstudiesareimportant toinvestigatetherobustnessofthemethoditselfwithrespecttosuch aspects.
5. Summaryandresearchperspectives
Weintroduceiterativemultiscalegradientcomputationalgorithms based on theiterativeMultiscale Finite Volume(i-MSFV)simulation method.By firstlyre-castingthei-MSFVin an algebraicframework, wederiveaflexible,multi-purposederivativecomputationframework whichaccountsfortheDirectandAdjointmethods.Theircomputational efficiencyisdiscussedanditisshownthattheyshareadvantages equiv-alenttothei-MSFVmethod.Thenewmethodsarevalidatedvia numer-icalexperimentsagainstnumericaldifferentiationareshowntoaddress thechallengesencounteredbytheMSFVgradientcomputation, specif-icallyassociatedwithhighheterogeneitycontrast.Fine-scale-accurate gradientsareobtainedforchallenginggeologicalmodelswithhigh per-meabilitycontrasts(e.g.theSPE-10comparativetestcase)ifthei-MSFV modelconvergestoanerrortoleranceof1.0𝑒−6.However,wehighlight
thatsuchaccuracyisnotnecessaryandonlyafewi-MSFV/smoothing it-erationsarerequiredforreasonablyaccurategradients.Furthercontrol ofthegradientqualityviaana-priorierrorestimatealongthe optimiza-tionprocessshouldallowforacomputationallyefficientoptimization method.Morespecifically,alessaccurategradientestimateusually suf-ficesintheearlyoptimizationiterations,whentheobjectivefunction convergenceismoresteep,whilemoreaccurategradientsarerequired intheflatterregionsforamoreprecisestoppingcriterion.Hence,the determinationofana-priorierrorestimateshouldbeakeydevelopment.
Acknowledgements
This work is part of Rafael Moraes’ PhD project, sponsored by
Petrobras’throughitsemployeedevelopmentprogram.Thediscussions heldwiththeDARSim(DelftAdvancedReservoirSimulation)members shallbeacknowledged.
References
Anterion, F., Eymard, R., Karcher, B., 1989. Use of parameter gradients for reservoir his- tory matching. SPE Symposium on Reservoir Simulation. Society of Petroleum Engi- neers. https://doi.org/10.2118/18433-MS .
Aziz, K. , et al. , 1993. Reservoir simulation grids: opportunities and problems. J. Petrol. Tech. 45 (07), 658–663 .
Cardoso, M. , Durlofsky, L. , Sarma, P. , 2009. Development and application of re- duced-order modeling procedures for subsurface flow simulation. Int. J. Numer. Method. Eng. 77 (9), 1322–1350 .
Chavent, G. , Dupuy, M. , Lemmonier, P. , 1975. History matching by use of optimal theory. Soc. Petrol. Eng. J. 15 (01), 74–86 .
Christie, M. , Blunt, M. , et al. , 2001. Tenth spe comparative solution project: A comparison of upscaling techniques. SPE Reservoir Simulation Symposium. Society of Petroleum Engineers .
van Doren, J.F., Markovinovi ć, R., Jansen, J.-D., 2006. Reduced-order optimal control of water flooding using proper orthogonal decomposition. Comp. Geosci. 10 (1), 137– 158. https://doi.org/10.1007/s10596-005-9014-2 .
Durlofsky, L.J. , 2005. Upscaling and gridding of fine scale geological models for flow simulation. 8th International Forum on Reservoir Simulation Iles Borromees, Stresa, Italy, 2024 .
Fonseca, R.M., Kahrobaei, S.S., Van Gastel, L.J.T., Leeuwenburgh, O., Jansen, J.D., 2015. Quantification of the impact of ensemble size on the quality of an ensemble gradient using principles of hypothesis testing. SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. https://doi.org/10.2118/173236-MS .
Frank, P., 2017. Iterative multiscale adjoint gradient computation for incompressible flow optimization in porous media.
Fu, J. , Caers, J. , Tchelepi, H.A. , 2011. A multiscale method for subsurface inverse model- ing: single-phase transient flow. Adv. Water Resour. 34 (8), 967–979 .
Fu, J. , Tchelepi, H.A. , Caers, J. , 2010. A multiscale adjoint method to compute sensitiv- ity coefficients for flow in heterogeneous porous media. Adv. Water Resour. 33 (6), 698–709 .
Hajibeygi, H. , 2011. Iterative multiscale finite volume method for multiphase flow in porous media with complex physics Ph.D. thesis. ETH Zurich .
Hajibeygi, H. , Bonfigli, G. , Hesse, M.A. , Jenny, P. , 2008. Iterative multiscale finite-volume method. J. Comput. Phys. 227 (19), 8604–8621 .
Hajibeygi, H. , Jenny, P. , 2011. Adaptive iterative multiscale finite volume method. J. Comput. Phys. 230 (3), 628–643 .
Hou, T.Y. , Wu, X.-H. , 1997. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys. 134 (1), 169–189 . Jansen, J.D. , 2011. Adjoint-based optimization of multi-phase flow through porous medi-
a–a review. Comp. Fluids 46 (1), 40–51 .
Jansen, J.D. , 2013. A simple algorithm to generate small geostatistical en- sembles for subsurface flow simulation. Research note.. Dept. of Geo- science and Engineering, Delft University of Technology, The Netherlands . uuid:6000459e-a0cb-40d1-843b-81650053e093.
Jansen, J. D., 2016. Gradient-based optimization of flow through porous media, v.3.
Jenny, P. , Lee, S.H. , Tchelepi, H.A. , 2003. Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. J. Comput. Phys. 187 (1), 47–67 .
Kraaijevanger, J.F.B.M., Egberts, P.J.P., Valstar, J.R., Buurman, H.W., 2007. Optimal wa- terflood design using the adjoint method. SPE Reservoir Simulation Symposium. So- ciety of Petroleum Engineers. https://doi.org/10.2118/105764-MS .
Krogstad, S. , Hauge, V.L. , Gulbransen, A. , 2011. Adjoint multiscale mixed finite elements. SPE J. 16 (01), 162–171 .
Li, R. , Reynolds, A. , Oliver, D. , et al. , 2003. History matching of three-phase flow produc- tion data. SPEJ 8 (04), 328–340 .
Lie, K.-A., Møyner, O., Natvig, J.R., Kozlova, A., Bratvedt, K., Watanabe, S., Li, Z., 2017. Successful application of multiscale methods in a real reservoir simulator en- vironment. Comp. Geosci. 21 (5), 981–998. https://doi.org/10.1007/s10596-017- 9627-2 .
Moraes, R. , Rodrigues, J. , Hajibeygi, H. , Jansen, J. , et al. , 2017. Multiscale gradient com- putation for multiphase flow in porous media. In: SPE Reservoir Simulation Confer- ence. Society of Petroleum Engineers .
de Moraes, R. J., Rodrigues, J. R., Hajibeygi, H., Jansen, J. D.,. Gradient computation for sequentially coupled subsurface optimization models. In Preparation.
de Moraes, R.J. , Rodrigues, J.R. , Hajibeygi, H. , Jansen, J.D. , 2017. Multiscale gradient computation for flow in heterogeneous porous media. J. Comput. Phys. 336, 644–663 . Oliver, D.S. , Reynolds, A.C. , Liu, N. , 2008. Inverse theory for petroleum reservoir charac-
terization and history matching. Cambridge University Press .
Remy, N. , Boucher, A. , Wu, J. , 2009. Applied Geostatistics with SGeMS: A User’s Guide. Cambridge University Press .
Rodrigues, J.R.P. , 2006. Calculating derivatives for automatic history matching. Comput. Geosci. 10 (1), 119–136 .
Saad, Y. , 2003. Iterative methods for sparse linear systems. SIAM .
Ţ ene, M. , Wang, Y. , Hajibeygi, H. , 2015. Adaptive algebraic multiscale solver for com- pressible flow in heterogeneous porous media. J. Comput. Phys. 300, 679–694 . Wang, Y. , Hajibeygi, H. , Tchelepi, H.A. , 2014. Algebraic multiscale solver for flow in
heterogeneous porous media. J. Comput. Phys. 259, 284–303 .
Zhou, H. , Tchelepi, H.A. , 2008. Operator-based multiscale method for compressible flow. SPE J. 13 (02), 523–539 .
