/[escript]/trunk/downunder/test/python/run_gravity.py
ViewVC logotype

Diff of /trunk/downunder/test/python/run_gravity.py

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

revision 4059 by caltinay, Wed Oct 31 01:01:24 2012 UTC revision 4060 by caltinay, Fri Nov 9 03:50:19 2012 UTC
# Line 22  __url__="https://launchpad.net/escript-f Line 22  __url__="https://launchpad.net/escript-f
22    
23  import logging  import logging
24  import os  import os
25    from esys.downunder import *
26  from esys.escript import unitsSI as U  from esys.escript import unitsSI as U
27  from esys.downunder.datasources import SyntheticDataSource,SmoothAnomaly  from esys.weipa import saveSilo
 from esys.downunder.inversions import GravityInversion  
28    
29    
30  try:  try:
# Line 32  try: Line 32  try:
32  except KeyError:  except KeyError:
33     WORKDIR='.'     WORKDIR='.'
34    
35  features=[SmoothAnomaly(lx=50*U.km, ly=20*U.km, lz=40*U.km, \  features=[SmoothAnomaly(lx=30*U.km, ly=20*U.km, lz=18.*U.km, \
36       x=100*U.km, y=3*U.km, depth=25*U.km, rho_inner=200., rho_outer=1e-6),\       x=8*U.km, y=3*U.km, depth=2.5*U.km, v_inner=200., v_outer=1e-6),\
37            SmoothAnomaly(lx=50*U.km, ly=20*U.km, lz=40*U.km,            SmoothAnomaly(lx=25*U.km, ly=20*U.km, lz=20*U.km,
38       x=400*U.km, y=1*U.km, depth=40*U.km, rho_inner=-200, rho_outer=1e-6)]       x=30*U.km, y=1*U.km, depth=18*U.km, v_inner=-200., v_outer=1e-6),\
39              SmoothAnomaly(lx=30*U.km, ly=20*U.km, lz=18.*U.km, \
40         x=68*U.km, y=3*U.km, depth=5*U.km, v_inner=200., v_outer=1e-6)]
41    
42  logger=logging.getLogger('inv')  logger=logging.getLogger('inv')
43  logger.setLevel(logging.INFO)  logger.setLevel(logging.INFO)
44  handler=logging.StreamHandler()  handler=logging.StreamHandler()
45  handler.setLevel(logging.INFO)  handler.setLevel(logging.INFO)
46  logger.addHandler(handler)  logger.addHandler(handler)
47  source=SyntheticDataSource(DIM=2, NE=20, l=500*U.km, h=60*U.km, features=features)  source=SyntheticData(DataSource.GRAVITY, DIM=2, NE=30, l=100*U.km, features=features)
48  source.setPadding(5)  domainbuilder=DomainBuilder(dim=2)
49    domainbuilder.addSource(source)
50    domainbuilder.setPadding(10)
51    domainbuilder.setVerticalExtents(depth=30*U.km, air_layer=10*U.km, num_cells=16)
52    
53  inv=GravityInversion()  inv=GravityInversion()
54  inv.setSolverTolerance(1e-5)  inv.setSolverTolerance(1e-9)
55  inv.setSolverMaxIterations(100)  inv.setSolverMaxIterations(100)
56  inv.setSolverOptions(initialHessian=100)  inv.setSolverOptions(initialHessian=10)
57  inv.setup(source)  inv.setWeights(mu_reg=1e-5)
58    inv.setup(domainbuilder)
59    
60  rho_new=inv.run()  rho_new=inv.run()
61    g, chi = inv.getForwardModel().getSurvey(0)
62    saveSilo(os.path.join(WORKDIR, 'gravinv'), density=rho_new, density_mask=inv.getRegularization().location_of_set_m, density_ref=source.getReferenceDensity(), g=g, chi=chi)
63    

Legend:
Removed from v.4059  
changed lines
  Added in v.4060

  ViewVC Help
Powered by ViewVC 1.1.26