/[escript]/trunk/finley/test/python/axisymm-splitB.py
ViewVC logotype

Diff of /trunk/finley/test/python/axisymm-splitB.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1638 by gross, Wed May 21 13:04:40 2008 UTC revision 1639 by gross, Mon Jul 14 08:55:25 2008 UTC
# Line 49  pressureStep2=LinearSinglePDE(dom) Line 49  pressureStep2=LinearSinglePDE(dom)
49  pressureStep2.setReducedOrderOn()  pressureStep2.setReducedOrderOn()
50  pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))  pressureStep2.setValue(q=whereZero(x[0]-L)+whereZero(x[1]-H))
51    
   
   
52  momentumStep3=LinearPDESystem(dom)  momentumStep3=LinearPDESystem(dom)
53  momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.])  momentumStep3.setValue(q=whereZero(x[0])*[1.,0.]+whereZero(x[1])*[0.,1.])
54  #  #
# Line 58  momentumStep3.setValue(q=whereZero(x[0]) Line 56  momentumStep3.setValue(q=whereZero(x[0])
56  #  #
57  U=Vector(0.,Solution(dom))  U=Vector(0.,Solution(dom))
58  p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3  p=ro*g*(L-ReducedSolution(dom).getX()[0])*(H-ReducedSolution(dom).getX()[1])/3
59  # p=ro*g*(H-ReducedSolution(dom).getX()[1])  p=ro*g*(H-ReducedSolution(dom).getX()[1])
60  stress=Tensor(0.,Function(dom))  dev_stress=Tensor(0.,Function(dom))
61    
62  t=dt  t=dt
63  istep=0  istep=0
# Line 78  while istep < nstep: Line 76  while istep < nstep:
76      sigma_d=(sign(inner(t_d,U))*alpha_w*t_d-n_d)*Pen*clip(inner(n_d,U),0.)      sigma_d=(sign(inner(t_d,U))*alpha_w*t_d-n_d)*Pen*clip(inner(n_d,U),0.)
77      print " sigma_d =",inf(sigma_d),sup(sigma_d)      print " sigma_d =",inf(sigma_d),sup(sigma_d)
78    
79      momentumStep1.setValue(D=ro*kronecker(dom),      momentumStep1.setValue(D=r*ro*kronecker(dom),
80                             Y=ro*U+dt*((stress[:,0]-p*teta3*kronecker(dom)[:,0])/r+[0.,-ro*g]),                             Y=r*ro*U+dt*r*[0.,-ro*g],
81                             # Y=r*ro*U+dt*r*[0.,-ro*g],                             X=-dt*r*(dev_stress-teta3*p*kronecker(dom)),
82                             X=-dt*(stress-teta3*p*kronecker(dom)),                             y=sigma_d*face_mask*r_b)
                            y=sigma_d*face_mask)  
83      U_star=momentumStep1.getSolution()      U_star=momentumStep1.getSolution()
84        saveVTK("u.xml",u=U_star,u0=U)
85        1/0
86      #      #
87      #  step 2:      #  step 2:
88      #      #
# Line 123  while istep < nstep: Line 122  while istep < nstep:
122      m=whereNegative(eta*gamma-alpha*p_pos)      m=whereNegative(eta*gamma-alpha*p_pos)
123      eta_d=m*eta+(1.-m)*alpha*p_pos/(gamma+small)        eta_d=m*eta+(1.-m)*alpha*p_pos/(gamma+small)  
124      print " viscosity =",inf(eta_d),sup(eta_d)      print " viscosity =",inf(eta_d),sup(eta_d)
125      stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))      dev_stress=eta_d*(symmetric(gg)-2./3.*vol*kronecker(dom))
126      #      #
127      # step size control:      # step size control:
128      #      #

Legend:
Removed from v.1638  
changed lines
  Added in v.1639

  ViewVC Help
Powered by ViewVC 1.1.26