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

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]
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)