/[escript]/trunk/downunder/py_src/datasources.py
ViewVC logotype

Diff of /trunk/downunder/py_src/datasources.py

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

revision 4115 by gross, Fri Dec 14 04:48:48 2012 UTC revision 4121 by gross, Wed Dec 19 00:24:50 2012 UTC
# Line 610  class SyntheticDataBase(DataSource): Line 610  class SyntheticDataBase(DataSource):
610      k=self.getReferenceProperty(domain)      k=self.getReferenceProperty(domain)
611      # calculate the corresponding potential      # calculate the corresponding potential
612      z=x[DIM-1]      z=x[DIM-1]
613      m_psi_ref=whereZero(z-sup(z))      m_psi_ref=whereZero(z-sup(z))+whereZero(z-inf(z))
614      if self.getDataType()==DataSource.GRAVITY:      if self.getDataType()==DataSource.GRAVITY:
615          pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)          pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)
616      else:      else:
617          pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)          pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)
618      pde.setSymmetryOn()      pde.setSymmetryOn()
619        #pde.getSolverOptions().setTolerance(1e-13)
620      psi_ref=pde.getSolution()      psi_ref=pde.getSolution()
621      del pde      del pde
622      if self.getDataType()==DataSource.GRAVITY:      if self.getDataType()==DataSource.GRAVITY:
623          data = -grad(psi_ref, ReducedFunction(domain))          data = -grad(psi_ref, ReducedFunction(domain))
624      else:      else:
625          data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))          data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))
626          
627      x=ReducedFunction(domain).getX()          x=ReducedFunction(domain).getX()    
628      if self.__full_knowledge:      if self.__full_knowledge:
629          sigma = whereNegative(x[DIM-1])          sigma = whereNegative(x[DIM-1])
# Line 712  class SyntheticData(SyntheticDataBase): Line 713  class SyntheticData(SyntheticDataBase):
713      """      """
714      defines synthetic  gravity/magnetic data based on harmonic property anomaly      defines synthetic  gravity/magnetic data based on harmonic property anomaly
715            
716          rho = mean + amplitude * sin(n_depth * pi /depth * z) * sin(n_length * pi /length * x - shift )          rho = mean + amplitude * sin(n_depth * pi /depth * (z+depth_offset)) * sin(n_length * pi /length * x - shift )
717                    
718      for all x and z<=0. for z>0 rho = 0.              for all x and z<=0. for z>0 rho = 0.        
719      """      """
720      def __init__(self, datatype,      def __init__(self, datatype,
721                         n_length=1,                         n_length=1,
722                         n_depth=1,                         n_depth=1,
723                         shift=0.,                         depth_offset=0.,
724                           depth=None,
725                         amplitude=None,                         amplitude=None,
726                         DIM=2,                         DIM=2,
727                         number_of_elements=10,                         number_of_elements=10,
# Line 734  class SyntheticData(SyntheticDataBase): Line 736  class SyntheticData(SyntheticDataBase):
736          :param n_length: number of oscillation in the anomaly data within the observation region.          :param n_length: number of oscillation in the anomaly data within the observation region.
737          :type n_length: ``int``          :type n_length: ``int``
738          :param n_depth: number of oscillation in the anomaly data below surface          :param n_depth: number of oscillation in the anomaly data below surface
739          :param shift: lateral phase shift in the  anomaly data          :param depth: vertical extend in the  anomaly data. If not present the depth of the domain is used.
740            :type depth: ``float``
741          :param amplitude: data amplitude. Default value is 200 *U.kg/U.m**3 for gravity and 0.1 for magnetic data.          :param amplitude: data amplitude. Default value is 200 *U.kg/U.m**3 for gravity and 0.1 for magnetic data.
742          :param features: list of features. It is recommended that the features do entirely lay below surface.          :param features: list of features. It is recommended that the features do entirely lay below surface.
743          :type features: ``list`` of ``SourceFeature``          :type features: ``list`` of ``SourceFeature``
# Line 761  class SyntheticData(SyntheticDataBase): Line 764  class SyntheticData(SyntheticDataBase):
764                                   spherical=spherical)                                   spherical=spherical)
765          self.__n_length = n_length          self.__n_length = n_length
766          self.__n_depth = n_depth          self.__n_depth = n_depth
767          self.__shift = shift          self.depth = depth
768            self.depth_offset=depth_offset
769          if amplitude == None:          if amplitude == None:
770          if datatype == DataSource.GRAVITY:          if datatype == DataSource.GRAVITY:
771              amplitude = 200 *U.kg/U.m**3              amplitude = 200 *U.kg/U.m**3
# Line 781  class SyntheticData(SyntheticDataBase): Line 785  class SyntheticData(SyntheticDataBase):
785              x=domain.getX()              x=domain.getX()
786              # set the reference data              # set the reference data
787              z=x[DIM-1]              z=x[DIM-1]
788              k=sin(self.__n_depth * np.pi/inf(z) * z) * whereNegative(z) * self.__amplitude              dd=self.depth
789                if  dd ==None :  dd=inf(z)
790                z2=(z+self.depth_offset)/(self.depth_offset-dd)
791                k=sin(self.__n_depth * np.pi  * z2) * whereNonNegative(z2) * whereNonPositive(z2-1.) * self.__amplitude
792              for i in xrange(DIM-1):              for i in xrange(DIM-1):
793             k*= sin(self.__n_length * np.pi /self.__length * x[i] - shift )             x_i=x[i]
794               m=whereNonNegative(x_i)*whereNonNegative(self.length-x_i)
795               k*= sin(self.__n_length*np.pi /self.length * x_i)*m
796              self._reference_data= k              self._reference_data= k
797          return self._reference_data          return self._reference_data
798    

Legend:
Removed from v.4115  
changed lines
  Added in v.4121

  ViewVC Help
Powered by ViewVC 1.1.26