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

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

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

revision 4853 by sshaw, Mon Mar 17 05:00:27 2014 UTC revision 4854 by gross, Wed Apr 9 06:17:09 2014 UTC
# Line 14  Line 14 
14  #  #
15  ##############################################################################  ##############################################################################
16    
17  """this a very simple test for inversion of acoustic data in the freqency domain"""  """this a very simple test for inversion of acoustic data in the freqency domain
18    
19    
20    
21    """
22    
23  __copyright__="""Copyright (c) 2003-2014 by University of Queensland  __copyright__="""Copyright (c) 2003-2014 by University of Queensland
24  http://www.uq.edu.au  http://www.uq.edu.au
# Line 42  logger.setLevel(logging.INFO) Line 46  logger.setLevel(logging.INFO)
46  # frequence range to be used for the inversion  # frequence range to be used for the inversion
47  OMEGA_MIN=1*U.Hz  OMEGA_MIN=1*U.Hz
48  OMEGA_MAX=10*U.Hz  OMEGA_MAX=10*U.Hz
49  N_OMEGA=3  N_OMEGA=1
50  # these are the taeget P-velocity and Q index  # these are the target P-velocity and Q index
51  V_P=2.*U.km  V_P=2.*U.km
52  Q=10  Q=10
53  # these are the reference values to be used  # these are the reference values to be used
# Line 67  domainbuilder.setVerticalExtents(depth=D Line 71  domainbuilder.setVerticalExtents(depth=D
71  # generate list of omegas:  # generate list of omegas:
72  omegas=np.linspace(OMEGA_MIN, OMEGA_MAX, num=N_OMEGA, endpoint=True)  omegas=np.linspace(OMEGA_MIN, OMEGA_MAX, num=N_OMEGA, endpoint=True)
73  # create data:  # create data:
74    #
75    #   - laplace u - o**2*SIGMA * u = 1.
76    #
77    #    -> u=1./o**2/SIGMA = data
78    #
79    #    to test the scaling an extra factor S is introduced.
80    #
81  SIGMA=1./( V_P * complex(1., -1./(2.*Q)))**2  SIGMA=1./( V_P * complex(1., -1./(2.*Q)))**2
82    SIGMA_prior=1./( V_P_prior * complex(1., -1./(2.*Q_prior)))**2
83  S={}  S={}
84  for o in omegas:  for o in omegas:
85      S[o]=complex(sqrt(o), sin(o))      S[o]=complex(sqrt(o)**2, sin(o))*0+1
86            
87      DD=-1./(SIGMA*o**2)*S[o]      DD=-1./(SIGMA*o**2)*S[o]
88      data=np.ones((NE_X+1,))*DD      data=np.ones((NE_X+1,))*DD
# Line 90  DIM=dom.getDim() Line 102  DIM=dom.getDim()
102        
103    # create regularization with two level set functions:    # create regularization with two level set functions:
104  z=Solution(dom).getX()[DIM-1]  z=Solution(dom).getX()[DIM-1]
105    x=Solution(dom).getX()[0]
106  reg_mask=Data(0.,(2,), Solution(dom))  reg_mask=Data(0.,(2,), Solution(dom))
107  reg_mask[0] = whereZero(z+inf(z))        # mask for locations where velocity and Q index are assumed to be known  reg_mask[0] = whereZero(z-inf(z))        # mask for locations where velocity and Q index are assumed to be known
108  reg_mask[1] = reg_mask[0]    reg_mask[1] = reg_mask[0]  
109    
110  regularization=Regularization(dom, numLevelSets=2,  regularization=Regularization(dom, numLevelSets=2,
# Line 129  for ds in domainbuilder.getTags(): Line 142  for ds in domainbuilder.getTags():
142      logger.debug("power for frequency %s Hz for source at %s added."%(ds.getFrequency, ds.getLocation()))      logger.debug("power for frequency %s Hz for source at %s added."%(ds.getFrequency, ds.getLocation()))
143      forwardmodels.append((acw, 0))      forwardmodels.append((acw, 0))
144    
145  m=AcousticVelocityMapping(V_P_prior, Q_prior)  vm=AcousticVelocityMapping(V_P, Q)
146    #
147    
148    z=Solution(dom).getX()[DIM-1]
149    x=Solution(dom).getX()[0]
150    z_bar=(z-inf(z))/(sup(z)-inf(z))
151    x_bar=(x-inf(x))/(sup(x)-inf(x))
152    s0=((SIGMA_prior.real-SIGMA.real)*z_bar+SIGMA.real)+SIGMA.real*0.1*sin((z_bar*x_bar)*15)
153    s1=((SIGMA_prior.imag-SIGMA.imag)*z_bar+SIGMA.imag)+SIGMA.imag*0.1*sin((z_bar*x_bar)*15)
154    m0=vm.getInverse(s0*[1,0]+s1*[0,1])
155    print "m0=",m0[0]
156    print "m1=",m0[1]
157    print "s0=",s0
158    print "s1=",s1
159    print "SIGMA=",SIGMA
160  # finally we can set up the cost function:  # finally we can set up the cost function:
161  cf=InversionCostFunction(regularization,  cf=InversionCostFunction(regularization,
162               (m,) , tuple(forwardmodels) )               (vm,) , tuple(forwardmodels) )
163    cf.setTradeOffFactorsRegularization(mu=[1e-1,1e-1], mu_c=None)
164                            
165  # and then start running the show:  # and then start running the show:
166  solver=MinimizerLBFGS()  solver=MinimizerLBFGS()
167  solver.setCostFunction(cf)  solver.setCostFunction(cf)
168  solver.setTolerance(1e-4)  solver.setTolerance(1e-6)
169  solver.setMaxIterations(50)  solver.setMaxIterations(180)
170  solver.run(Data(0.,(2,),Solution(dom)))  solver.run(m0)
171  m=solver.getResult()  m=solver.getResult()
172  sigma = cf.getProperties(m)  sigma = cf.getProperties(m)
173        
174  print sigma  print "m0=",m[0]
175    print "m1=",m[1]
176    print "s0=",sigma[0]
177    print "s1=",sigma[1]
178    print "SIGMA=",SIGMA
179    
 1/0  
 #=============================================================  
 # interesting parameters:  
 n_humps_h = 3  
 n_humps_v = 1  
 mu = 100  
 n_cells_in_data = 100  
 # ignore:  
 full_knowledge = False  
 depth_offset = 0. * U.km  
 #  
 DIM = 2  
 n_cells_in_data = max(n_humps_h*7, n_cells_in_data)  
 l_data = 100. * U.km  
 l_pad = 40. * U.km  
 THICKNESS = 20. * U.km  
 l_air = 20. * U.km  
 n_cells_v = max(  
         int((2*l_air+THICKNESS+depth_offset)/l_data*n_cells_in_data + 0.5), 25)  
   
   
 source=SyntheticData(DataSource.GRAVITY, n_length=n_humps_h, n_depth=n_humps_v,  
         depth=THICKNESS+depth_offset, depth_offset=depth_offset,  
         DIM=DIM, number_of_elements=n_cells_in_data, length=l_data,  
         data_offset=0, full_knowledge=full_knowledge)  
   
 domainbuilder=DomainBuilder(dim=DIM)  
 domainbuilder.addSource(source)  
 domainbuilder.setVerticalExtents(depth=l_air+THICKNESS+depth_offset,  
                                  air_layer=l_air, num_cells=n_cells_v)  
 domainbuilder.setPadding(l_pad)  
 domainbuilder.fixDensityBelow(depth=THICKNESS+depth_offset)  
   
 inv=GravityInversion()  
 inv.setSolverTolerance(1e-4)  
 inv.setSolverMaxIterations(50)  
 inv.setup(domainbuilder)  
 inv.getCostFunction().setTradeOffFactorsModels(mu)  
   
   
 rho_new = inv.run()  
 rho_ref = source.getReferenceProperty()  
 print("rho_new = %s"%rho_new)  
 print("rho = %s"%rho_ref)  
 g, chi = inv.getCostFunction().getForwardModel().getSurvey(0)  
 saveSilo(os.path.join(WORKDIR, 'results_gravity_2d'), density=rho_new, density_ref=rho_ref, g=g, chi=chi)  
180    

Legend:
Removed from v.4853  
changed lines
  Added in v.4854

  ViewVC Help
Powered by ViewVC 1.1.26