48 |
import util |
import util |
49 |
import math |
import math |
50 |
|
|
51 |
|
##### Added by Artak |
52 |
|
from Numeric import zeros,Int,Float32,Float64 |
53 |
|
################################### |
54 |
|
|
55 |
|
|
56 |
class TimeIntegrationManager: |
class TimeIntegrationManager: |
57 |
""" |
""" |
58 |
a simple mechanism to manage time dependend values. |
a simple mechanism to manage time dependend values. |
544 |
|
|
545 |
return x,r |
return x,r |
546 |
|
|
547 |
|
|
548 |
|
############################ |
549 |
|
# Added by Artak |
550 |
|
#################################3 |
551 |
|
|
552 |
|
#Apply a sequence of k Givens rotations, used within gmres codes |
553 |
|
# vrot=givapp(c, s, vin, k) |
554 |
|
def givapp(c,s,vin,k): |
555 |
|
vrot=vin |
556 |
|
for i in range(k+1): |
557 |
|
w1=c[i]*vrot[i]-s[i]*vrot[i+1] |
558 |
|
w2=s[i]*vrot[i]+c[i]*vrot[i+1] |
559 |
|
vrot[i:i+2]=w1,w2 |
560 |
|
return vrot |
561 |
|
|
562 |
|
def GMRES(b, Aprod, Msolve, bilinearform, stoppingcriterium, x=None, iter_max=100): |
563 |
|
|
564 |
|
from numarray import dot |
565 |
|
|
566 |
|
v0=b |
567 |
|
iter=0 |
568 |
|
if x==None: |
569 |
|
x=0*b |
570 |
|
else: |
571 |
|
b += (-1.)*Aprod(x) |
572 |
|
r=b |
573 |
|
|
574 |
|
rhat=Msolve(r) |
575 |
|
|
576 |
|
rhat_dot_r = bilinearform(rhat, rhat) |
577 |
|
norm_r=math.sqrt(rhat_dot_r) |
578 |
|
|
579 |
|
if rhat_dot_r<0: raise NegativeNorm,"negative norm." |
580 |
|
|
581 |
|
h=zeros((iter_max,iter_max),Float32) |
582 |
|
c=zeros(iter_max,Float32) |
583 |
|
s=zeros(iter_max,Float32) |
584 |
|
g=zeros(iter_max,Float32) |
585 |
|
v=[] |
586 |
|
|
587 |
|
v.append(rhat/norm_r) |
588 |
|
rho=norm_r |
589 |
|
g[0]=rho |
590 |
|
|
591 |
|
while not stoppingcriterium(rho,rho,norm_r): |
592 |
|
|
593 |
|
if iter >= iter_max: raise MaxIterReached,"maximum number of %s steps reached."%iter_max |
594 |
|
|
595 |
|
|
596 |
|
vhat=Aprod(v[iter]) |
597 |
|
p=Msolve(vhat) |
598 |
|
|
599 |
|
v.append(p) |
600 |
|
|
601 |
|
v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1])) |
602 |
|
|
603 |
|
# Modified Gram-Schmidt |
604 |
|
for j in range(iter+1): |
605 |
|
h[j][iter]=bilinearform(v[j],v[iter+1]) |
606 |
|
v[iter+1]+=(-1.)*h[j][iter]*v[j] |
607 |
|
|
608 |
|
h[iter+1][iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1])) |
609 |
|
v_norm2=h[iter+1][iter] |
610 |
|
|
611 |
|
|
612 |
|
# Reorthogonalize if needed |
613 |
|
if v_norm1 + 0.001*v_norm2 == v_norm1: #Brown/Hindmarsh condition (default) |
614 |
|
for j in range(iter+1): |
615 |
|
hr=bilinearform(v[j],v[iter+1]) |
616 |
|
h[j][iter]=h[j][iter]+hr #vhat |
617 |
|
v[iter+1] +=(-1.)*hr*v[j] |
618 |
|
|
619 |
|
v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1])) |
620 |
|
h[iter+1][iter]=v_norm2 |
621 |
|
|
622 |
|
# watch out for happy breakdown |
623 |
|
if v_norm2 != 0: |
624 |
|
v[iter+1]=v[iter+1]/h[iter+1][iter] |
625 |
|
|
626 |
|
# Form and store the information for the new Givens rotation |
627 |
|
if iter > 0 : |
628 |
|
hhat=[] |
629 |
|
for i in range(iter+1) : hhat.append(h[i][iter]) |
630 |
|
hhat=givapp(c[0:iter],s[0:iter],hhat,iter-1); |
631 |
|
for i in range(iter+1) : h[i][iter]=hhat[i] |
632 |
|
|
633 |
|
mu=math.sqrt(h[iter][iter]*h[iter][iter]+h[iter+1][iter]*h[iter+1][iter]) |
634 |
|
if mu!=0 : |
635 |
|
c[iter]=h[iter][iter]/mu |
636 |
|
s[iter]=-h[iter+1][iter]/mu |
637 |
|
h[iter][iter]=c[iter]*h[iter][iter]-s[iter]*h[iter+1][iter] |
638 |
|
h[iter+1][iter]=0.0 |
639 |
|
g[iter:iter+2]=givapp(c[iter],s[iter],g[iter:iter+2],0) |
640 |
|
|
641 |
|
# Update the residual norm |
642 |
|
rho=abs(g[iter+1]) |
643 |
|
iter+=1 |
644 |
|
|
645 |
|
# At this point either iter > iter_max or rho < tol. |
646 |
|
# It's time to compute x and leave. |
647 |
|
|
648 |
|
if iter > 0 : |
649 |
|
y=zeros(iter,Float32) |
650 |
|
y[iter-1] = g[iter-1] / h[iter-1][iter-1] |
651 |
|
if iter > 1 : |
652 |
|
i=iter-2 |
653 |
|
while i>=0 : |
654 |
|
y[i] = ( g[i] - dot(h[i][i+1:iter], y[i+1:iter])) / h[i][i] |
655 |
|
i=i-1 |
656 |
|
xhat=v[iter-1]*y[iter-1] |
657 |
|
for i in range(iter-1): |
658 |
|
xhat += v[i]*y[i] |
659 |
|
else : xhat=v[0] |
660 |
|
|
661 |
|
x += xhat |
662 |
|
|
663 |
|
return x |
664 |
|
|
665 |
|
############################################# |
666 |
|
|
667 |
class ArithmeticTuple(object): |
class ArithmeticTuple(object): |
668 |
""" |
""" |
669 |
tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float. |
tuple supporting inplace update x+=y and scaling x=a*y where x,y is an ArithmeticTuple and a is a float. |
733 |
out.append(other*self[i]) |
out.append(other*self[i]) |
734 |
return ArithmeticTuple(*tuple(out)) |
return ArithmeticTuple(*tuple(out)) |
735 |
|
|
736 |
|
######################### |
737 |
|
# Added by Artak |
738 |
|
######################### |
739 |
|
def __div__(self,other): |
740 |
|
""" |
741 |
|
dividing from the right |
742 |
|
|
743 |
|
@param other: scaling factor |
744 |
|
@type other: C{float} |
745 |
|
@return: itemwise self/other |
746 |
|
@rtype: L{ArithmeticTuple} |
747 |
|
""" |
748 |
|
out=[] |
749 |
|
for i in range(len(self)): |
750 |
|
out.append(self[i]/other) |
751 |
|
return ArithmeticTuple(*tuple(out)) |
752 |
|
|
753 |
|
def __rdiv__(self,other): |
754 |
|
""" |
755 |
|
dividing from the left |
756 |
|
|
757 |
|
@param other: scaling factor |
758 |
|
@type other: C{float} |
759 |
|
@return: itemwise other/self |
760 |
|
@rtype: L{ArithmeticTuple} |
761 |
|
""" |
762 |
|
out=[] |
763 |
|
for i in range(len(self)): |
764 |
|
out.append(other/self[i]) |
765 |
|
return ArithmeticTuple(*tuple(out)) |
766 |
|
|
767 |
|
##########################################33 |
768 |
|
|
769 |
def __iadd__(self,other): |
def __iadd__(self,other): |
770 |
""" |
""" |
771 |
in-place add of other to self |
in-place add of other to self |
848 |
def __inner(self,p,r): |
def __inner(self,p,r): |
849 |
return self.inner(p,r[1]) |
return self.inner(p,r[1]) |
850 |
|
|
851 |
|
def __inner_p(self,p1,p2): |
852 |
|
return self.inner(p1,p2) |
853 |
|
|
854 |
def __stoppingcriterium(self,norm_r,r,p): |
def __stoppingcriterium(self,norm_r,r,p): |
855 |
return self.stoppingcriterium(r[1],r[0],p) |
return self.stoppingcriterium(r[1],r[0],p) |
856 |
|
|
857 |
|
def __stoppingcriterium_GMRES(self,norm_r,rho,r): |
858 |
|
return self.stoppingcriterium_GMRES(rho,r) |
859 |
|
|
860 |
def setTolerance(self,tolerance=1.e-8): |
def setTolerance(self,tolerance=1.e-8): |
861 |
self.__tol=tolerance |
self.__tol=tolerance |
862 |
def getTolerance(self): |
def getTolerance(self): |
866 |
def getSubProblemTolerance(self): |
def getSubProblemTolerance(self): |
867 |
return self.__reduction*self.getTolerance() |
return self.__reduction*self.getTolerance() |
868 |
|
|
869 |
def solve(self,v,p,max_iter=20, verbose=False, show_details=False): |
def solve(self,v,p,max_iter=20, verbose=False, show_details=False, solver='GMRES'): |
870 |
""" |
""" |
871 |
solves the saddle point problem using initial guesses v and p. |
solves the saddle point problem using initial guesses v and p. |
872 |
|
|
877 |
|
|
878 |
# Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask) |
# Az=f is solved as A(z-v)=f-Av (z-v = 0 on fixed_u_mask) |
879 |
|
|
880 |
|
|
881 |
self.__z=v+self.solve_A(v,p*0) |
self.__z=v+self.solve_A(v,p*0) |
882 |
|
|
883 |
Bz=self.B(self.__z) |
Bz=self.B(self.__z) |
884 |
# |
# |
885 |
# solve BA^-1B^*p = Bz |
# solve BA^-1B^*p = Bz |
890 |
# A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask) |
# A(v-z)=Az-B^*p-Az = f -Az - B^*p (v-z=0 on fixed_u_mask) |
891 |
# |
# |
892 |
self.iter=0 |
self.iter=0 |
893 |
if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter |
if solver=='GMRES': |
894 |
p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p) |
if self.verbose: print "enter GMRES method (iter_max=%s)"%max_iter |
895 |
return r[0],p |
p=GMRES(Bz,self.__Aprod_GMRES,self.__Msolve_GMRES,self.__inner_p,self.__stoppingcriterium_GMRES,iter_max=max_iter, x=p*1) |
896 |
|
u=v+self.solve_A(v,p) |
897 |
|
|
898 |
|
else: |
899 |
|
if self.verbose: print "enter PCG method (iter_max=%s)"%max_iter |
900 |
|
p,r=PCG(ArithmeticTuple(self.__z*1.,Bz),self.__Aprod,self.__Msolve,self.__inner,self.__stoppingcriterium,iter_max=max_iter, x=p*1) |
901 |
|
u=r[0] |
902 |
|
u=v+self.solve_A(v,p) # Lutz to check !!!!! |
903 |
|
|
904 |
|
return u,p |
905 |
|
|
906 |
def __Msolve(self,r): |
def __Msolve(self,r): |
907 |
return self.solve_prec(r[1]) |
return self.solve_prec(r[1]) |
908 |
|
|
909 |
|
def __Msolve_GMRES(self,r): |
910 |
|
return self.solve_prec(r) |
911 |
|
|
912 |
|
|
913 |
def __Aprod(self,p): |
def __Aprod(self,p): |
914 |
# return BA^-1B*p |
# return BA^-1B*p |
915 |
#solve Av =-B^*p as Av =f-Az-B^*p |
#solve Av =-B^*p as Av =f-Az-B^*p |
916 |
v=self.solve_A(self.__z,p) |
v=self.solve_A(self.__z,p) |
917 |
return ArithmeticTuple(v, self.B(v)) |
return ArithmeticTuple(v, self.B(v)) |
918 |
|
|
919 |
|
def __Aprod_GMRES(self,p): |
920 |
|
# return BA^-1B*p |
921 |
|
#solve Av =-B^*p as Av =f-Az-B^*p |
922 |
|
v=self.solve_A(self.__z,p) |
923 |
|
return self.B(v) |
924 |
|
|
925 |
class SaddlePointProblem(object): |
class SaddlePointProblem(object): |
926 |
""" |
""" |
1177 |
return util.whereNonZero(util.interpolate(out,function_space)) |
return util.whereNonZero(util.interpolate(out,function_space)) |
1178 |
|
|
1179 |
|
|
1180 |
|
|