# Diff of /trunk/escript/py_src/flows.py

revision 1517 by artak, Fri Apr 18 02:36:37 2008 UTC revision 2445 by gross, Fri May 29 03:23:25 2009 UTC
# Line 1  Line 1
1  # \$Id:\$  ########################################################
2  #  #
3  #######################################################  # Copyright (c) 2003-2008 by University of Queensland
4    # Earth Systems Science Computational Center (ESSCC)
5    # http://www.uq.edu.au/esscc
6  #  #
7  #       Copyright 2008 by University of Queensland  # Primary Business: Queensland, Australia
#
#######################################################
10  #  #
11    ########################################################
12
14    Earth Systems Science Computational Center (ESSCC)
15    http://www.uq.edu.au/esscc
20
21  """  """
22  Some models for flow  Some models for flow
# Line 24  Some models for flow Line 30  Some models for flow
30  """  """
31
32  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
http://www.access.edu.au
__url__="http://www.iservo.edu.au/esys"
__version__="\$Revision:\$"
__date__="\$Date:\$"
33
34  from escript import *  from escript import *
35  import util  import util
36  from linearPDEs import LinearPDE  from linearPDEs import LinearPDE, LinearPDESystem, LinearSinglePDE
37  from pdetools import HomogeneousSaddlePointProblem  from pdetools import HomogeneousSaddlePointProblem,Projector, ArithmeticTuple, PCG, NegativeNorm, GMRES
38
39    class DarcyFlow(object):
40        """
41        solves the problem
42
43        M{u_i+k_{ij}*p_{,j} = g_i}
44        M{u_{i,i} = f}
45
46        where M{p} represents the pressure and M{u} the Darcy flux. M{k} represents the permeability,
47
48        @note: The problem is solved in a least squares formulation.
49        """
50
51        def __init__(self, domain, weight=None, useReduced=False):
52            """
53            initializes the Darcy flux problem
54            @param domain: domain of the problem
55            @type domain: L{Domain}
56            """
57            self.domain=domain
58            if weight == None:
59               s=self.domain.getSize()
60               self.__l=(3.*util.longestEdge(self.domain)*s/util.sup(s))**2
61            else:
62               self.__l=weight
63            self.__pde_v=LinearPDESystem(domain)
64            if useReduced: self.__pde_v.setReducedOrderOn()
65            self.__pde_v.setSymmetryOn()
66            self.__pde_v.setValue(D=util.kronecker(domain), A=self.__l*util.outer(util.kronecker(domain),util.kronecker(domain)))
67            # self.__pde_v.setSolverMethod(preconditioner=self.__pde_v.ILU0)
68            self.__pde_p=LinearSinglePDE(domain)
69            self.__pde_p.setSymmetryOn()
70            if useReduced: self.__pde_p.setReducedOrderOn()
71            self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
72            self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
73            self.setTolerance()
74            self.setAbsoluteTolerance()
75            self.setSubProblemTolerance()
76
77        def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None):
78            """
79            assigns values to model parameters
80
81            @param f: volumetic sources/sinks
82            @type f: scalar value on the domain (e.g. L{Data})
83            @param g: flux sources/sinks
84            @type g: vector values on the domain (e.g. L{Data})
85            @param location_of_fixed_pressure: mask for locations where pressure is fixed
86            @type location_of_fixed_pressure: scalar value on the domain (e.g. L{Data})
87            @param location_of_fixed_flux:  mask for locations where flux is fixed.
88            @type location_of_fixed_flux: vector values on the domain (e.g. L{Data})
89            @param permeability: permeability tensor. If scalar C{s} is given the tensor with
90                                 C{s} on the main diagonal is used. If vector C{v} is given the tensor with
91                                 C{v} on the main diagonal is used.
92            @type permeability: scalar, vector or tensor values on the domain (e.g. L{Data})
93
94            @note: the values of parameters which are not set by calling C{setValue} are not altered.
95            @note: at any point on the boundary of the domain the pressure (C{location_of_fixed_pressure} >0)
96                   or the normal component of the flux (C{location_of_fixed_flux[i]>0} if direction of the normal
97                   is along the M{x_i} axis.
98            """
99            if f !=None:
100               f=util.interpolate(f, self.__pde_v.getFunctionSpaceForCoefficient("X"))
101               if f.isEmpty():
102                   f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X"))
103               else:
104                   if f.getRank()>0: raise ValueError,"illegal rank of f."
105               self.__f=f
106            if g !=None:
107               g=util.interpolate(g, self.__pde_p.getFunctionSpaceForCoefficient("Y"))
108               if g.isEmpty():
109                 g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y"))
110               else:
111                 if not g.getShape()==(self.domain.getDim(),):
112                   raise ValueError,"illegal shape of g"
113               self.__g=g
114
115            if location_of_fixed_pressure!=None: self.__pde_p.setValue(q=location_of_fixed_pressure)
116            if location_of_fixed_flux!=None: self.__pde_v.setValue(q=location_of_fixed_flux)
117
118            if permeability!=None:
119               perm=util.interpolate(permeability,self.__pde_p.getFunctionSpaceForCoefficient("A"))
120               if perm.getRank()==0:
121                   perm=perm*util.kronecker(self.domain.getDim())
122               elif perm.getRank()==1:
123                   perm, perm2=Tensor(0.,self.__pde_p.getFunctionSpaceForCoefficient("A")), perm
124                   for i in range(self.domain.getDim()): perm[i,i]=perm2[i]
125               elif perm.getRank()==2:
126                  pass
127               else:
128                  raise ValueError,"illegal rank of permeability."
129               self.__permeability=perm
130               self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability))
131
132        def setTolerance(self,rtol=1e-4):
133            """
134            sets the relative tolerance C{rtol} used to terminate the solution process. The iteration is terminated if
135
136            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
137
138            where C{atol} is an absolut tolerance (see L{setAbsoluteTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
139
140            @param rtol: relative tolerance for the pressure
141            @type rtol: non-negative C{float}
142            """
143            if rtol<0:
144                raise ValueError,"Relative tolerance needs to be non-negative."
145            self.__rtol=rtol
146        def getTolerance(self):
147            """
148            returns the relative tolerance
149
150            @return: current relative tolerance
151            @rtype: C{float}
152            """
153            return self.__rtol
154
155        def setAbsoluteTolerance(self,atol=0.):
156            """
157            sets the absolute tolerance C{atol} used to terminate the solution process. The iteration is terminated if
158
159            M{|g-v-Qp| <= atol + rtol * min( max( |g-v|, |Qp| ), max( |v|, |g-Qp| ) ) }
160
161            where C{rtol} is an absolut tolerance (see L{setTolerance}), M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
162
163            @param atol: absolute tolerance for the pressure
164            @type atol: non-negative C{float}
165            """
166            if atol<0:
167                raise ValueError,"Absolute tolerance needs to be non-negative."
168            self.__atol=atol
169        def getAbsoluteTolerance(self):
170           """
171           returns the absolute tolerance
172
173           @return: current absolute tolerance
174           @rtype: C{float}
175           """
176           return self.__atol
177
178        def setSubProblemTolerance(self,rtol=None):
179             """
180             Sets the relative tolerance to solve the subproblem(s). If C{rtol} is not present
181             C{self.getTolerance()**2} is used.
182
183             @param rtol: relative tolerence
184             @type rtol: positive C{float}
185             """
186             if rtol == None:
187                  if self.getTolerance()<=0.:
188                      raise ValueError,"A positive relative tolerance must be set."
189                  self.__sub_tol=max(util.EPSILON**(0.75),self.getTolerance()**2)
190             else:
191                 if rtol<=0:
192                     raise ValueError,"sub-problem tolerance must be positive."
193                 self.__sub_tol=max(util.EPSILON**(0.75),rtol)
194
195        def getSubProblemTolerance(self):
196             """
197             Returns the subproblem reduction factor.
198
199             @return: subproblem reduction factor
200             @rtype: C{float}
201             """
202             return self.__sub_tol
203
204        def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, max_num_corrections=10):
205             """
206             solves the problem.
207
208             The iteration is terminated if the residual norm is less then self.getTolerance().
209
210             @param u0: initial guess for the flux. At locations in the domain marked by C{location_of_fixed_flux} the value of C{u0} is kept unchanged.
211             @type u0: vector value on the domain (e.g. L{Data}).
212             @param p0: initial guess for the pressure. At locations in the domain marked by C{location_of_fixed_pressure} the value of C{p0} is kept unchanged.
213             @type p0: scalar value on the domain (e.g. L{Data}).
214             @param verbose: if set some information on iteration progress are printed
215             @type verbose: C{bool}
216             @param show_details:  if set information on the subiteration process are printed.
217             @type show_details: C{bool}
218             @return: flux and pressure
219             @rtype: C{tuple} of L{Data}.
220
221             @note: The problem is solved as a least squares form
222
223             M{(I+D^*D)u+Qp=D^*f+g}
224             M{Q^*u+Q^*Qp=Q^*g}
225
226             where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}.
227             We eliminate the flux form the problem by setting
228
229             M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} with u=u0 on location_of_fixed_flux
230
231             form the first equation. Inserted into the second equation we get
232
233             M{Q^*(I-(I+D^*D)^{-1})Qp= Q^*(g-(I+D^*D)^{-1}(D^*f+g))} with p=p0  on location_of_fixed_pressure
234
235             which is solved using the PCG method (precondition is M{Q^*Q}). In each iteration step
236             PDEs with operator M{I+D^*D} and with M{Q^*Q} needs to be solved using a sub iteration scheme.
237             """
238             self.verbose=verbose
239             self.show_details= show_details and self.verbose
240             rtol=self.getTolerance()
241             atol=self.getAbsoluteTolerance()
242             if self.verbose: print "DarcyFlux: initial sub tolerance = %e"%self.getSubProblemTolerance()
243
244             num_corrections=0
245             converged=False
246             p=p0
247             norm_r=None
248             while not converged:
249                   v=self.getFlux(p, fixed_flux=u0, show_details=self.show_details)
250                   Qp=self.__Q(p)
251                   norm_v=self.__L2(v)
252                   norm_Qp=self.__L2(Qp)
253                   if norm_v == 0.:
254                      if norm_Qp == 0.:
255                         return v,p
256                      else:
257                        fac=norm_Qp
258                   else:
259                      if norm_Qp == 0.:
260                        fac=norm_v
261                      else:
262                        fac=2./(1./norm_v+1./norm_Qp)
263                   ATOL=(atol+rtol*fac)
264                   if self.verbose:
265                        print "DarcyFlux: L2 norm of v = %e."%norm_v
266                        print "DarcyFlux: L2 norm of k.grad(p) = %e."%norm_Qp
267                        print "DarcyFlux: absolute tolerance ATOL = %e."%ATOL
268                   if norm_r == None or norm_r>ATOL:
269                       if num_corrections>max_num_corrections:
270                             raise ValueError,"maximum number of correction steps reached."
271                       p,r, norm_r=PCG(self.__g-util.interpolate(v,Function(self.domain))-Qp,self.__Aprod,p,self.__Msolve_PCG,self.__inner_PCG,atol=0.5*ATOL, rtol=0.,iter_max=max_iter, verbose=self.verbose)
272                       num_corrections+=1
273                   else:
274                       converged=True
275             return v,p
276        def __L2(self,v):
277             return util.sqrt(util.integrate(util.length(util.interpolate(v,Function(self.domain)))**2))
278
279        def __Q(self,p):
281
282        def __Aprod(self,dp):
283              self.__pde_v.setTolerance(self.getSubProblemTolerance())
284              if self.show_details: print "DarcyFlux: Applying operator"
285              Qdp=self.__Q(dp)
286              self.__pde_v.setValue(Y=-Qdp,X=Data(), r=Data())
287              du=self.__pde_v.getSolution(verbose=self.show_details, iter_max = 100000)
288              # self.__pde_v.getOperator().saveMM("proj.mm")
289              return Qdp+du
290        def __inner_GMRES(self,r,s):
291             return util.integrate(util.inner(r,s))
292
293        def __inner_PCG(self,p,r):
294             return util.integrate(util.inner(self.__Q(p), r))
295
296        def __Msolve_PCG(self,r):
297              self.__pde_p.setTolerance(self.getSubProblemTolerance())
298              if self.show_details: print "DarcyFlux: Applying preconditioner"
299              self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r), Y=Data(), r=Data())
300              # self.__pde_p.getOperator().saveMM("prec.mm")
301              return self.__pde_p.getSolution(verbose=self.show_details, iter_max = 100000)
302
303        def getFlux(self,p=None, fixed_flux=Data(), show_details=False):
304            """
305            returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux}
306            on locations where C{location_of_fixed_flux} is positive (see L{setValue}).
307            Note that C{g} and C{f} are used, see L{setValue}.
308
309            @param p: pressure.
310            @type p: scalar value on the domain (e.g. L{Data}).
311            @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}.
312            @type fixed_flux: vector values on the domain (e.g. L{Data}).
313            @param tol: relative tolerance to be used.
314            @type tol: positive C{float}.
315            @return: flux
316            @rtype: L{Data}
317            @note: the method uses the least squares solution M{u=(I+D^*D)^{-1}(D^*f-g-Qp)} where M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}}
318                   for the permeability M{k_{ij}}
319            """
320            self.__pde_v.setTolerance(self.getSubProblemTolerance())
321            g=self.__g
322            f=self.__f
323            self.__pde_v.setValue(X=self.__l*f*util.kronecker(self.domain), r=fixed_flux)
324            if p == None:
325               self.__pde_v.setValue(Y=g)
326            else:
327               self.__pde_v.setValue(Y=g-self.__Q(p))
328            return self.__pde_v.getSolution(verbose=show_details, iter_max=100000)
329
331        """       """
332        solves       solves
333
334            -(eta*(u_{i,j}+u_{j,i}))_j - p_i = f_i            -(eta*(u_{i,j}+u_{j,i}))_j + p_i = f_i-stress_{ij,j}
335                  u_{i,i}=0                  u_{i,i}=0
336
338            eta*(u_{i,j}+u_{j,i})*n_j=surface_stress            eta*(u_{i,j}+u_{j,i})*n_j-p*n_i=surface_stress +stress_{ij}n_j
339
340        if surface_stress is not give 0 is assumed.       if surface_stress is not given 0 is assumed.
341
342        typical usage:       typical usage:
343
344              sp=StokesProblemCartesian(domain)              sp=StokesProblemCartesian(domain)
345              sp.setTolerance()              sp.setTolerance()
346              sp.initialize(...)              sp.initialize(...)
347              v,p=sp.solve(v0,p0)              v,p=sp.solve(v0,p0)
348        """       """
349        def __init__(self,domain,**kwargs):       def __init__(self,domain,**kwargs):
350             """
351             initialize the Stokes Problem
352
353             @param domain: domain of the problem. The approximation order needs to be two.
354             @type domain: L{Domain}
355             @warning: The apprximation order needs to be two otherwise you may see oscilations in the pressure.
356             """
358           self.domain=domain           self.domain=domain
359           self.vol=util.integrate(1.,Function(self.domain))           self.vol=util.integrate(1.,Function(self.domain))
360           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())           self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
361           self.__pde_u.setSymmetryOn()           self.__pde_u.setSymmetryOn()
362  #         self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)           # self.__pde_u.setSolverMethod(self.__pde_u.DIRECT)
363                         # self.__pde_u.setSolverMethod(preconditioner=LinearPDE.ILU0)
364
365           self.__pde_prec=LinearPDE(domain)           self.__pde_prec=LinearPDE(domain)
366           self.__pde_prec.setReducedOrderOn()           self.__pde_prec.setReducedOrderOn()
367             # self.__pde_prec.setSolverMethod(self.__pde_prec.LUMPING)
368           self.__pde_prec.setSymmetryOn()           self.__pde_prec.setSymmetryOn()
369
370           self.__pde_proj=LinearPDE(domain)           self.__pde_proj=LinearPDE(domain)
371           self.__pde_proj.setReducedOrderOn()           self.__pde_proj.setReducedOrderOn()
372         self.__pde_proj.setValue(D=1)
373           self.__pde_proj.setSymmetryOn()           self.__pde_proj.setSymmetryOn()
self.__pde_proj.setValue(D=1.)
374
377            """
378            assigns values to the model parameters
379
380            @param f: external force
381            @type f: L{Vector} object in L{FunctionSpace} L{Function} or similar
383            @type fixed_u_mask: L{Vector} object on L{FunctionSpace} L{Solution} or similar
384            @param eta: viscosity
385            @type eta: L{Scalar} object on L{FunctionSpace} L{Function} or similar
386            @param surface_stress: normal surface stress
387            @type eta: L{Vector} object on L{FunctionSpace} L{FunctionOnBoundary} or similar
388            @param stress: initial stress
389        @type stress: L{Tensor} object on L{FunctionSpace} L{Function} or similar
390            @note: All values needs to be set.
391
392            """
393          self.eta=eta          self.eta=eta
394          A =self.__pde_u.createCoefficientOfGeneralPDE("A")          A =self.__pde_u.createCoefficient("A")
395      self.__pde_u.setValue(A=Data())      self.__pde_u.setValue(A=Data())
396          for i in range(self.domain.getDim()):          for i in range(self.domain.getDim()):
397          for j in range(self.domain.getDim()):          for j in range(self.domain.getDim()):
398              A[i,j,j,i] += 1.              A[i,j,j,i] += 1.
399              A[i,j,i,j] += 1.              A[i,j,i,j] += 1.
400      self.__pde_prec.setValue(D=1./self.eta)      self.__pde_prec.setValue(D=1/self.eta)
402            self.__f=f
403            self.__surface_stress=surface_stress
404            self.__stress=stress
405
406         def Bv(self,v):
407             """
408             returns inner product of element p and div(v)
409
410             @param p: a pressure increment
411             @param v: a residual
412             @return: inner product of element p and div(v)
413             @rtype: C{float}
414             """
415             self.__pde_proj.setValue(Y=-util.div(v))
416             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
417             return self.__pde_proj.getSolution()
418
419         def inner_pBv(self,p,Bv):
420             """
421             returns inner product of element p and Bv=-div(v)
422
423             @param p: a pressure increment
424             @param v: a residual
425             @return: inner product of element p and Bv=-div(v)
426             @rtype: C{float}
427             """
428             return util.integrate(util.interpolate(p,Function(self.domain))*util.interpolate(Bv,Function(self.domain)))
429
430         def inner_p(self,p0,p1):
431             """
432             Returns inner product of p0 and p1
433
434        def B(self,arg):           @param p0: a pressure
435           d=util.div(arg)           @param p1: a pressure
436           self.__pde_proj.setValue(Y=d)           @return: inner product of p0 and p1
437           self.__pde_proj.setTolerance(self.getSubProblemTolerance())           @rtype: C{float}
438           return self.__pde_proj.getSolution(verbose=self.show_details)           """
439             s0=util.interpolate(p0/self.eta,Function(self.domain))
440        def inner(self,p0,p1):           s1=util.interpolate(p1/self.eta,Function(self.domain))
s0=util.interpolate(p0,Function(self.domain))
s1=util.interpolate(p1,Function(self.domain))
441           return util.integrate(s0*s1)           return util.integrate(s0*s1)
442
443        def getStress(self,u):       def norm_v(self,v):
445           return 2.*self.eta*util.symmetric(mg)           returns the norm of v
446
447        def solve_A(self,u,p):           @param v: a velovity
448             @return: norm of v
449             @rtype: non-negative C{float}
450           """           """
452
453         def getV(self, p, v0):
454             """
455             return the value for v for a given p (overwrite)
456
457             @param p: a pressure
458             @param v0: a initial guess for the value v to return.
459             @return: v given as M{v= A^{-1} (f-B^*p)}
460           """           """
461           self.__pde_u.setTolerance(self.getSubProblemTolerance())           self.__pde_u.setTolerance(self.getSubProblemTolerance())
462           self.__pde_u.setValue(X=-self.getStress(u)-p*util.kronecker(self.domain))           self.__pde_u.setValue(Y=self.__f, y=self.__surface_stress, r=v0)
463           return  self.__pde_u.getSolution(verbose=self.show_details)           if self.__stress.isEmpty():
464                self.__pde_u.setValue(X=p*util.kronecker(self.domain))
465             else:
466                self.__pde_u.setValue(X=self.__stress+p*util.kronecker(self.domain))
467             out=self.__pde_u.getSolution(verbose=self.show_details)
468             return  out
469
470         def norm_Bv(self,Bv):
471            """
472            Returns Bv (overwrite).
473
474            @rtype: equal to the type of p
475            @note: boundary conditions on p should be zero!
476            """
477            return util.sqrt(util.integrate(util.interpolate(Bv,Function(self.domain))**2))
478
479        def solve_prec(self,p):       def solve_AinvBt(self,p):
480           self.__pde_prec.setTolerance(self.getSubProblemTolerance())           """
481           self.__pde_prec.setValue(Y=p)           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}
q=self.__pde_prec.getSolution(verbose=self.show_details)
return q
def stoppingcriterium(self,Bv,v,p):
n_r=util.sqrt(self.inner(Bv,Bv))
n_v=util.Lsup(v)
if self.verbose: print "PCG step %s: L2(div(v)) = %s, Lsup(v)=%s"%(self.iter,n_r,n_v)
self.iter+=1
if n_r <= self.vol**(1./2.-1./self.domain.getDim())*n_v*self.getTolerance():
if self.verbose: print "PCG terminated after %s steps."%self.iter
return True
else:
return False
def stoppingcriterium2(self,norm_r,norm_b,solver='GMRES'):
if self.verbose: print "%s step %s: L2(r) = %s, L2(b)*TOL=%s"%(solver,self.iter,norm_r,norm_b*self.getTolerance())
self.iter+=1
if norm_r <= norm_b*self.getTolerance():
if self.verbose: print "%s terminated after %s steps."%(solver,self.iter)
return True
else:
return False
482
483             @param p: a pressure increment
484             @return: the solution of M{Av=B^*p}
485             @note: boundary conditions on v should be zero!
486             """
487             self.__pde_u.setTolerance(self.getSubProblemTolerance())
488             self.__pde_u.setValue(Y=Data(), y=Data(), r=Data(),X=-p*util.kronecker(self.domain))
489             out=self.__pde_u.getSolution(verbose=self.show_details)
490             return  out
491
492         def solve_prec(self,Bv):
493             """
494             applies preconditioner for for M{BA^{-1}B^*} to M{Bv}
495             with accuracy L{self.getSubProblemTolerance()}
496
497             @param v: velocity increment
498             @return: M{p=P(Bv)} where M{P^{-1}} is an approximation of M{BA^{-1}B^*}
499             @note: boundary conditions on p are zero.
500             """
501             self.__pde_prec.setValue(Y=Bv)
502             self.__pde_prec.setTolerance(self.getSubProblemTolerance())
503             return self.__pde_prec.getSolution(verbose=self.show_details)
504
505    # rename solve_prec and change argument v to Bv
506    # chnage the argument of inner_pBv to v->Bv