/[escript]/branches/symbolic_from_3470/ripley/test/python/gravity_gamma.py
ViewVC logotype

Diff of /branches/symbolic_from_3470/ripley/test/python/gravity_gamma.py

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

revision 3876 by caltinay, Sun Mar 18 23:30:57 2012 UTC revision 3877 by gross, Mon Mar 19 00:52:00 2012 UTC
# Line 59  else: Line 59  else:
59     x_cord=domain.getX()-[L0, L1, H_earth]     x_cord=domain.getX()-[L0, L1, H_earth]
60  LL=max(L,H)  LL=max(L,H)
61        
62  m_psi_ref=whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) +  whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))  m_psi=whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) +  whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))
63  for i in range(DIM-1):  for i in range(DIM-1):
64      m_psi_ref= m_psi_ref + whereZero(x_cord[i]-inf(x_cord[i])) +  whereZero(x_cord[i]-sup(x_cord[i]))      m_psi= m_psi + whereZero(x_cord[i]-inf(x_cord[i])) +  whereZero(x_cord[i]-sup(x_cord[i]))
65  #m_rho=wherePositive(domain.getX()[DIM-1]-H_earth)  #m_rho=wherePositive(domain.getX()[DIM-1]-H_earth)
66  #m_rho=whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) +  whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))  #m_rho=whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) +  whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))
67  m_rho=wherePositive(domain.getX()[DIM-1]-H_earth) + whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) +  whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))  m_rho=wherePositive(domain.getX()[DIM-1]-H_earth) + whereZero(x_cord[DIM-1]-inf(x_cord[DIM-1])) +  whereZero(x_cord[DIM-1]-sup(x_cord[DIM-1]))
# Line 78  for f in feastures: Line 78  for f in feastures:
78    
79  # get the reference potential:  # get the reference potential:
80  pde=LinearSinglePDE(domain)  pde=LinearSinglePDE(domain)
81  pde.setValue(A=kronecker(domain), Y=4*pi*rho_ref, q=m_psi_ref)  pde.setValue(A=kronecker(domain), Y=4*pi*rho_ref, q=m_psi)
82  pde.getSolverOptions().setVerbosityOn()  pde.getSolverOptions().setVerbosityOn()
83  pde.setSymmetryOn()  pde.setSymmetryOn()
84  #pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)  #pde.getSolverOptions().setSolverMethod(pde.getSolverOptions().DIRECT)
# Line 87  del pde Line 87  del pde
87  d_obs=kronecker(DIM)[DIM-1]  d_obs=kronecker(DIM)[DIM-1]
88  g_hat=grad(psi_ref)[DIM-1]  g_hat=grad(psi_ref)[DIM-1]
89  beta=1/1000000.  beta=1/1000000.
90  beta=1/100.  #beta=1/100.
91  #  #
92  #  where do we know the gravity:  #  where do we know the gravity:
93  #  #
# Line 109  print "====== Use variational problem  = Line 109  print "====== Use variational problem  =
109  psi_s=Symbol("psi", (), dim=DIM)  psi_s=Symbol("psi", (), dim=DIM)
110  rho_s=Symbol("rho", (), dim=DIM)  rho_s=Symbol("rho", (), dim=DIM)
111  gamma_s=Symbol("gamma", (), dim=DIM)  gamma_s=Symbol("gamma", (), dim=DIM)
112    gamma_s=0.5
113    gamma_s=0.6
114  #g=Symbol("g", (), dim=DIM)  #g=Symbol("g", (), dim=DIM)
115    
116  v=VariationalProblem(domain, u=psi_s,p=rho_s, debug=VariationalProblem.DEBUG3)  v=VariationalProblem(domain, u=psi_s,p=rho_s, debug=VariationalProblem.DEBUG3)
117  v.setValue( H = 0.5*chi*(grad(psi_s)[DIM-1]-g_hat)**2 + beta/gamma_s * (length(grad(rho_s))**2 + EPSILON)**gamma_s,  v.setValue( H = 0.5*chi*(grad(psi_s)[DIM-1]-g_hat)**2 + beta/gamma_s * (length(grad(rho_s))**2 + EPSILON**2)**gamma_s,
118              X=grad(psi_s), Y=-4*pi*rho_s/LL**2,              X=grad(psi_s), Y=-4*pi*rho_s/LL**2,
119              qp=m_rho, q=m_psi_ref)              qp=m_rho, q=m_psi)
120  v.getNonlinearPDE().getLinearSolverOptions().setSolverMethod(v.getNonlinearPDE().getLinearSolverOptions().DIRECT)  v.getNonlinearPDE().getLinearSolverOptions().setSolverMethod(v.getNonlinearPDE().getLinearSolverOptions().DIRECT)
121                            
122  rho_v, psi_v, lag=v.getSolution(psi=0, rho=1, gamma=1./2)  # gamma=1 is the interesting case!  rho_v, psi_v, lag=v.getSolution(psi=0, rho=1, gamma=2)  # gamma=0.5 is the interesting case!
123  print "rho =",rho_v  print "rho =",rho_v
124  print "rho_ref =",rho_ref*g0  print "rho_ref =",rho_ref*g0
125  print "psi =",psi_v  print "psi =",psi_v

Legend:
Removed from v.3876  
changed lines
  Added in v.3877

  ViewVC Help
Powered by ViewVC 1.1.26