--- trunk/escript/py_src/flows.py 2008/12/15 05:09:02 2156 +++ trunk/escript/py_src/flows.py 2009/01/12 06:37:07 2208 @@ -48,7 +48,7 @@ @note: The problem is solved in a least squares formulation. """ - def __init__(self, domain): + def __init__(self, domain,useReduced=False): """ initializes the Darcy flux problem @param domain: domain of the problem @@ -56,12 +56,15 @@ """ self.domain=domain self.__pde_v=LinearPDESystem(domain) - self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain))) + if useReduced: self.__pde_v.setReducedOrderOn() self.__pde_v.setSymmetryOn() + self.__pde_v.setValue(D=util.kronecker(domain), A=util.outer(util.kronecker(domain),util.kronecker(domain))) self.__pde_p=LinearSinglePDE(domain) self.__pde_p.setSymmetryOn() + if useReduced: self.__pde_p.setReducedOrderOn() self.__f=Scalar(0,self.__pde_v.getFunctionSpaceForCoefficient("X")) self.__g=Vector(0,self.__pde_v.getFunctionSpaceForCoefficient("Y")) + self.__ATOL= None def setValue(self,f=None, g=None, location_of_fixed_pressure=None, location_of_fixed_flux=None, permeability=None): """ @@ -119,42 +122,125 @@ self.__pde_p.setValue(A=util.transposed_tensor_mult(self.__permeability,self.__permeability)) - def getFlux(self,p, fixed_flux=Data(),tol=1.e-8, show_details=False): + def getFlux(self,p=None, fixed_flux=Data(),tol=1.e-8, show_details=False): """ returns the flux for a given pressure C{p} where the flux is equal to C{fixed_flux} on locations where C{location_of_fixed_flux} is positive (see L{setValue}). - Note that C{g} and C{f} are used, L{setValue}. + Note that C{g} and C{f} are used, see L{setValue}. @param p: pressure. @type p: scalar value on the domain (e.g. L{Data}). @param fixed_flux: flux on the locations of the domain marked be C{location_of_fixed_flux}. @type fixed_flux: vector values on the domain (e.g. L{Data}). @param tol: relative tolerance to be used. - @type tol: positive float. + @type tol: positive C{float}. @return: flux @rtype: L{Data} @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}} for the permeability M{k_{ij}} """ self.__pde_v.setTolerance(tol) - self.__pde_v.setValue(Y=self.__g, X=self.__f*util.kronecker(self.domain), r=fixed_flux) + g=self.__g + f=self.__f + self.__pde_v.setValue(X=f*util.kronecker(self.domain), r=fixed_flux) + if p == None: + self.__pde_v.setValue(Y=g) + else: + self.__pde_v.setValue(Y=g-util.tensor_mult(self.__permeability,util.grad(p))) return self.__pde_v.getSolution(verbose=show_details) - def solve(self,u0,p0,atol=0,rtol=1e-8, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8): + def getPressure(self,v=None, fixed_pressure=Data(),tol=1.e-8, show_details=False): + """ + returns the pressure for a given flux C{v} where the pressure is equal to C{fixed_pressure} + on locations where C{location_of_fixed_pressure} is positive (see L{setValue}). + Note that C{g} is used, see L{setValue}. + + @param v: flux. + @type v: vector-valued on the domain (e.g. L{Data}). + @param fixed_pressure: pressure on the locations of the domain marked be C{location_of_fixed_pressure}. + @type fixed_pressure: vector values on the domain (e.g. L{Data}). + @param tol: relative tolerance to be used. + @type tol: positive C{float}. + @return: pressure + @rtype: L{Data} + @note: the method uses the least squares solution M{p=(Q^*Q)^{-1}Q^*(g-u)} where and M{(Qp)_i=k_{ij}p_{,j}} + for the permeability M{k_{ij}} + """ + self.__pde_v.setTolerance(tol) + g=self.__g + self.__pde_p.setValue(r=fixed_pressure) + if v == None: + self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g-v)) + else: + self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,g)) + return self.__pde_p.getSolution(verbose=show_details) + + def setTolerance(self,atol=0,rtol=1e-8,p_ref=None,v_ref=None): + """ + set the tolerance C{ATOL} used to terminate the solution process. It is used + + M{ATOL = atol + rtol * max( |g-v_ref|, |Qp_ref| )} + + where M{|f|^2 = integrate(length(f)^2)} and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}}. If C{v_ref} or C{p_ref} is not present zero is assumed. + + The iteration is terminated if for the current approximation C{p}, flux C{v=(I+D^*D)^{-1}(D^*f-g-Qp)} and their residual + + M{r=Q^*(g-Qp-v)} + + the condition + + M{<(Q^*Q)^{-1} r,r> <= ATOL} + + holds. M{D} is the M{div} operator and M{(Qp)_i=k_{ij}p_{,j}} for the permeability M{k_{ij}} + + @param atol: absolute tolerance for the pressure + @type atol: non-negative C{float} + @param rtol: relative tolerance for the pressure + @type rtol: non-negative C{float} + @param p_ref: reference pressure. If not present zero is used. You may use physical arguments to set a resonable value for C{p_ref}, use the + L{getPressure} method or use the value from a previous time step. + @type p_ref: scalar value on the domain (e.g. L{Data}). + @param v_ref: reference velocity. If not present zero is used. You may use physical arguments to set a resonable value for C{v_ref}, use the + L{getFlux} method or use the value from a previous time step. + @type v_ref: vector-valued on the domain (e.g. L{Data}). + @return: used absolute tolerance. + @rtype: positive C{float} + """ + g=self.__g + if not v_ref == None: + f1=util.integrate(util.length(util.interpolate(g-v_ref,Function(self.domain)))**2) + else: + f1=util.integrate(util.length(util.interpolate(g))**2) + if not p_ref == None: + f2=util.integrate(util.length(util.tensor_mult(self.__permeability,util.grad(p_ref)))**2) + else: + f2=0 + self.__ATOL= atol + rtol * util.sqrt(max(f1,f2)) + if self.__ATOL<=0: + raise ValueError,"Positive tolerance (=%e) is expected."%self.__ATOL + return self.__ATOL + + def getTolerance(self): + """ + returns the current tolerance. + + @return: used absolute tolerance. + @rtype: positive C{float} + """ + if self.__ATOL==None: + raise ValueError,"no tolerance is defined." + return self.__ATOL + + def solve(self,u0,p0, max_iter=100, verbose=False, show_details=False, sub_rtol=1.e-8): """ solves the problem. - The iteration is terminated if the error in the pressure is less then C{rtol * |q| + atol} where - C{|q|} denotes the norm of the right hand side (see escript user's guide for details). + The iteration is terminated if the residual norm is less then self.getTolerance(). @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. @type u0: vector value on the domain (e.g. L{Data}). @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. @type p0: scalar value on the domain (e.g. L{Data}). - @param atol: absolute tolerance for the pressure - @type atol: non-negative C{float} - @param rtol: relative tolerance for the pressure - @type rtol: non-negative C{float} @param sub_rtol: tolerance to be used in the sub iteration. It is recommended that M{sub_rtol",f0-f1 + return f0-f1 def __Msolve_PCG(self,r): if self.show_details: print "DarcyFlux: Applying preconditioner" - self.__pde_p.setValue(X=-util.transposed_tensor_mult(self.__permeability,r[0]+r[1])) + self.__pde_p.setValue(X=util.transposed_tensor_mult(self.__permeability,r[0]-r[1]), r=Data()) return self.__pde_p.getSolution(verbose=self.show_details) class StokesProblemCartesian(HomogeneousSaddlePointProblem):