/[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 4125 by gross, Wed Jan 2 06:15:00 2013 UTC revision 4132 by caltinay, Fri Jan 11 05:46:49 2013 UTC
# Line 516  class SmoothAnomaly(SourceFeature): Line 516  class SmoothAnomaly(SourceFeature):
516    
517  ##############################################################################  ##############################################################################
518  class SyntheticDataBase(DataSource):  class SyntheticDataBase(DataSource):
519    """      """
520    Base class to define reference data based on a given property distribution (density or      Base class to define reference data based on a given property distribution
521    susceptibility). Data are collected from a square region of vertical extend `length` on a      (density or susceptibility). Data are collected from a square region of
522    grid with ``number_of_elements`` cells in each direction.      vertical extent `length` on a grid with `number_of_elements` cells in each
523          direction.
524    The synthetic data are constructed by solving the appropriate forward problem. Data can be sampled      
525    with an offset from the surface at z=0 or using the entire subsurface region.      The synthetic data are constructed by solving the appropriate forward
526    """      problem. Data can be sampled with an offset from the surface at z=0 or
527    def __init__(self, datatype,      using the entire subsurface region.
528          DIM=2,      """
529          number_of_elements=10,      def __init__(self, datatype,
530          length=1*U.km,                   DIM=2,
531          B_b=None,                   number_of_elements=10,
532          data_offset=0,                   length=1*U.km,
533          full_knowledge=False,                   B_b=None,
534          spherical=False):                   data_offset=0,
535      """                   full_knowledge=False,
536                     spherical=False):
537            """
538          :param datatype: data type indicator          :param datatype: data type indicator
539          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
540          :param DIM: spatial dimension          :param DIM: number of spatial dimension
541          :type DIM: 2 or 3          :type DIM: ``int`` (2 or 3)
542          :param number_of_elements: lateral number of elements in the region where data are collected          :param number_of_elements: lateral number of elements in the region
543                                       where data are collected
544          :type number_of_elements: ``int``          :type number_of_elements: ``int``
545          :param length: lateral extend of the region where data are collected          :param length: lateral extend of the region where data are collected
546          :type length: ``float``          :type length: ``float``
# Line 545  class SyntheticDataBase(DataSource): Line 548  class SyntheticDataBase(DataSource):
548          :type B_b: ``list`` of ``Scalar``          :type B_b: ``list`` of ``Scalar``
549          :param data_offset: offset of the data collection region from the surface          :param data_offset: offset of the data collection region from the surface
550          :type data_offset: ``float``          :type data_offset: ``float``
551          :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.          :param full_knowledge: if ``True`` data are collected from the entire
552                                   subsurface region. This is mainly for testing.
553          :type full_knowledge: ``Bool``          :type full_knowledge: ``Bool``
554          :param spherical: if ``True`` sperical coordinates are used (ignored)          :param spherical: if ``True`` sperical coordinates are used (ignored for now)
555          :type spherical: ``Bool``          :type spherical: ``Bool``
556      """          """
557      super(SyntheticDataBase,self).__init__()          super(SyntheticDataBase,self).__init__()
558      if not datatype in [self.GRAVITY,self.MAGNETIC]:          if not datatype in [self.GRAVITY,self.MAGNETIC]:
559          raise ValueError("Invalid value for datatype parameter")                    raise ValueError("Invalid value for datatype parameter")      
560      self.DIM=DIM          self.DIM=DIM
561      self.number_of_elements=number_of_elements          self.number_of_elements=number_of_elements
562      self.length=length          self.length=length
563      self.__datatype = datatype          self.__datatype = datatype
564                
565      self.__spherical = spherical            self.__spherical = spherical  
566      self.__full_knowledge= full_knowledge          self.__full_knowledge= full_knowledge
567      self.__data_offset=data_offset          self.__data_offset=data_offset
568      self.__B_b =None          self.__B_b =None
569      # this is for Cartesian (FIXME ?)          # this is for Cartesian (FIXME ?)
570      if datatype  ==  self.MAGNETIC:          if datatype  ==  self.MAGNETIC:
571          if self.DIM<3:              if self.DIM < 3:
572            self.__B_b =  np.array([-B_b[2],  -B_b[0]])                  self.__B_b = np.array([-B_b[2], -B_b[0]])
573          else:              else:
574            self.__B_b = ([-B_b[1],  -B_b[2],  -B_b[0]])                      self.__B_b = ([-B_b[1], -B_b[2], -B_b[0]])    
575      self.__origin = [0]*(DIM-1)          self.__origin = [0]*(DIM-1)
576      self.__delta = [float(length)/number_of_elements]*(DIM-1)          self.__delta = [float(length)/number_of_elements]*(DIM-1)
577      self.__nPts = [number_of_elements]*(DIM-1)          self.__nPts = [number_of_elements]*(DIM-1)
578      self._reference_data=None          self._reference_data=None
579                
580    def getDataExtents(self):      def getDataExtents(self):
581      """          """
582      returns the lateral data extend of the data set          returns the lateral data extend of the data set
583      """          """
584      return (list(self.__origin), list(self.__nPts), list(self.__delta))          return (list(self.__origin), list(self.__nPts), list(self.__delta))
585    
586    def getDataType(self):      def getDataType(self):
587      """          """
588      returns the data type          returns the data type
589      """          """
590      return self.__datatype          return self.__datatype
591                
592    def getSurveyData(self, domain, origin, number_of_elements, spacing):      def getSurveyData(self, domain, origin, number_of_elements, spacing):
593      """          """
594      returns the survey data placed on a given domain.          returns the survey data placed on a given domain.
595                
596      :param domain: domain on which the data are to be placed          :param domain: domain on which the data are to be placed
597      :type param: ``Domain``          :type domain: ``Domain``
598      :param origin: origin of the domain          :param origin: origin of the domain
599      :type origin: ``list`` of ``float``          :type origin: ``list`` of ``float``
600      :param number_of_elements: number of elements (or cells) in each spatial direction used          :param number_of_elements: number of elements (or cells) in each
601                    span the domain                                     spatial direction used to span the domain
602      :type number_of_elements: `list`` of ``int``          :type number_of_elements: ``list`` of ``int``
603      :param spacing: cell size in each spatial direction          :param spacing: cell size in each spatial direction
604      :type spacing: ``list`` of ``float``          :type spacing: ``list`` of ``float``
605      :return: observed gravity field or magnetic flux density for each cell in the domain and          :return: observed gravity field or magnetic flux density for each cell
606      for each cell an indicator 1/0 if the data are valid or not.                   in the domain and for each cell an indicator 1/0 if the data
607      :rtype: pair of ``Scalar``                   are valid or not.
608      """          :rtype: pair of ``Scalar``
609      pde=LinearSinglePDE(domain)          """
610      DIM=domain.getDim()          pde=LinearSinglePDE(domain)
611      x=domain.getX()          DIM=domain.getDim()
612      # set the reference data          x=domain.getX()
613                # set the reference data
614      k=self.getReferenceProperty(domain)          
615      # calculate the corresponding potential          k=self.getReferenceProperty(domain)
616      z=x[DIM-1]          # calculate the corresponding potential
617      m_psi_ref=whereZero(z-sup(z))          z=x[DIM-1]
618      if self.getDataType()==DataSource.GRAVITY:          m_psi_ref=whereZero(z-sup(z))
619          pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)          if self.getDataType()==DataSource.GRAVITY:
620      else:              pde.setValue(A=kronecker(domain), Y=-4*np.pi*U.Gravitational_Constant*self._reference_data, q=m_psi_ref)
621          pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)          else:
622      pde.setSymmetryOn()              pde.setValue(A=kronecker(domain), X=self._reference_data*self.__B_b, q=m_psi_ref)
623      #pde.getSolverOptions().setTolerance(1e-13)          pde.setSymmetryOn()
624      psi_ref=pde.getSolution()          #pde.getSolverOptions().setTolerance(1e-13)
625      del pde          psi_ref=pde.getSolution()
626      if self.getDataType()==DataSource.GRAVITY:          del pde
627          data = -grad(psi_ref, ReducedFunction(domain))          if self.getDataType()==DataSource.GRAVITY:
628      else:              data = -grad(psi_ref, ReducedFunction(domain))
629          data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))          else:
630                data = self._reference_data*self.__B_b-grad(psi_ref, ReducedFunction(domain))
631    
632      x=ReducedFunction(domain).getX()              x=ReducedFunction(domain).getX()    
633      if self.__full_knowledge:          if self.__full_knowledge:
634          sigma = whereNegative(x[DIM-1])              sigma = whereNegative(x[DIM-1])
635      else:          else:
636                      sigma=1.
637          sigma=1.              # limit mask to non-padding in horizontal area        
638          # limit mask to non-padding in horizontal area                      for i in range(DIM-1):
639          for i in range(DIM-1):                  x_i=x[i]
640          x_i=x[i]                  sigma=sigma * wherePositive(x_i) * whereNegative(x_i-(sup(x_i)+inf(x_i)))
641          sigma=sigma * wherePositive(x_i) * whereNegative(x_i-(sup(x_i)+inf(x_i)))              # limit mask to one cell thickness at z=0
642          # limit mask to one cell thickness at z=0              z=x[DIM-1]
643          z=x[DIM-1]              oo=int(self.__data_offset/spacing[DIM-1]+0.5)*spacing[DIM-1]
644          oo=int(self.__data_offset/spacing[DIM-1]+0.5)*spacing[DIM-1]              sigma = sigma * whereNonNegative(z-oo) * whereNonPositive(z-oo-spacing[DIM-1])
645          sigma = sigma * whereNonNegative(z-oo) * whereNonPositive(z-oo-spacing[DIM-1])          return data,sigma
646      return data,sigma          
647            def getReferenceProperty(self, domain=None):
648    def getReferenceProperty(self, domain=None):          """
649      """          Returns the reference density Data object that was used to generate
650      Returns the reference density Data object that was used to generate          the gravity/susceptibility anomaly data.
651      the gravity/susceptibility anomaly data.          
652                :return: the density or susceptibility anomaly used to create the
653      :return: the density or susceptibility anomaly used to create the survey data.                   survey data
654      :note: it can be assumed that in the first call the ``domain`` argument is present so the          :note: it can be assumed that in the first call the ``domain``
655      actual anomaly data can be created. In subsequent calls this may not be true.                 argument is present so the actual anomaly data can be created.
656      :note: method needs to be overwritten                 In subsequent calls this may not be true.
657      """          :note: method needs to be overwritten
658      raise NotImplementedError()                """
659            raise NotImplementedError()      
660                
661  class SyntheticFeatureData(SyntheticDataBase):  class SyntheticFeatureData(SyntheticDataBase):
662      """      """
# Line 696  class SyntheticFeatureData(SyntheticData Line 702  class SyntheticFeatureData(SyntheticData
702    
703      def getReferenceProperty(self, domain=None):      def getReferenceProperty(self, domain=None):
704          """          """
705      Returns the reference density Data object that was used to generate          Returns the reference density Data object that was used to generate
706      the gravity/susceptibility anomaly data.          the gravity/susceptibility anomaly data.
707          """          """
708          if self._reference_data == None:          if self._reference_data == None:
709              DIM=domain.getDim()              DIM=domain.getDim()
# Line 733  class SyntheticData(SyntheticDataBase): Line 739  class SyntheticData(SyntheticDataBase):
739          """          """
740          :param datatype: data type indicator          :param datatype: data type indicator
741          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``          :type datatype: ``DataSource.GRAVITY``, ``DataSource.MAGNETIC``
742          :param n_length: number of oscillation in the anomaly data within the observation region.          :param n_length: number of oscillations in the anomaly data within the
743                             observation region
744          :type n_length: ``int``          :type n_length: ``int``
745          :param n_depth: number of oscillation in the anomaly data below surface          :param n_depth: number of oscillations in the anomaly data below surface
746          :param depth: vertical extend in the  anomaly data. If not present the depth of the domain is used.          :param depth: vertical extent in the anomaly data. If not present the
747                          depth of the domain is used.
748          :type depth: ``float``          :type depth: ``float``
749          :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
750          :param features: list of features. It is recommended that the features do entirely lay below surface.                            gravity and 0.1 for magnetic data.
751          :type features: ``list`` of ``SourceFeature``          :param DIM: spatial dimensionality
752          :param DIM: spatial dimension          :type DIM: ``int`` (2 or 3)
753          :type DIM: 2 or 3          :param number_of_elements: lateral number of elements in the region
754          :param number_of_elements: lateral number of elements in the region where data are collected                                     where data are collected
755          :type number_of_elements: ``int``          :type number_of_elements: ``int``
756          :param length: lateral extend of the region where data are collected          :param length: lateral extent of the region where data are collected
757          :type length: ``float``          :type length: ``float``
758          :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude]. Only used for magnetic data.          :param B_b: background magnetic flux density [B_r, B_latiude, B_longitude].
759                        Only used for magnetic data.
760          :type B_b: ``list`` of ``Scalar``          :type B_b: ``list`` of ``Scalar``
761          :param data_offset: offset of the data collection region from the surface          :param data_offset: offset of the data collection region from the surface
762          :type data_offset: ``float``          :type data_offset: ``float``
763          :param full_knowledge: if ``True`` data are collected from the entire subsurface region. This is mainly for testing.          :param full_knowledge: if ``True`` data are collected from the entire
764                                   subsurface region. This is mainly for testing.
765          :type full_knowledge: ``Bool``          :type full_knowledge: ``Bool``
766          :param spherical: if ``True`` sperical coordinates are used (ignored)          :param spherical: if ``True`` sperical coordinates are used (ignored for now)
767          :type spherical: ``Bool``          :type spherical: ``Bool``
768          """                """      
769          super(SyntheticData,self).__init__(          super(SyntheticData,self).__init__(
770                                   datatype=datatype, DIM=DIM, number_of_elements=number_of_elements,                                   datatype=datatype, DIM=DIM,
771                                     number_of_elements=number_of_elements,
772                                   length=length, B_b=B_b,                                   length=length, B_b=B_b,
773                                   data_offset=data_offset,                                   data_offset=data_offset,
774                                   full_knowledge=full_knowledge,                                   full_knowledge=full_knowledge,
# Line 767  class SyntheticData(SyntheticDataBase): Line 778  class SyntheticData(SyntheticDataBase):
778          self.depth = depth          self.depth = depth
779          self.depth_offset=depth_offset          self.depth_offset=depth_offset
780          if amplitude == None:          if amplitude == None:
781          if datatype == DataSource.GRAVITY:              if datatype == DataSource.GRAVITY:
782              amplitude = 200 *U.kg/U.m**3                  amplitude = 200 *U.kg/U.m**3
783          else:              else:
784               amplitude =0.1                  amplitude = 0.1
785          self.__amplitude = amplitude          self.__amplitude = amplitude
786    
   
   
787      def getReferenceProperty(self, domain=None):      def getReferenceProperty(self, domain=None):
788          """          """
789          Returns the reference density Data object that was used to generate          Returns the reference density Data object that was used to generate
# Line 790  class SyntheticData(SyntheticDataBase): Line 799  class SyntheticData(SyntheticDataBase):
799              z2=(z+self.depth_offset)/(self.depth_offset-dd)              z2=(z+self.depth_offset)/(self.depth_offset-dd)
800              k=sin(self.__n_depth * np.pi  * z2) * whereNonNegative(z2) * whereNonPositive(z2-1.) * self.__amplitude              k=sin(self.__n_depth * np.pi  * z2) * whereNonNegative(z2) * whereNonPositive(z2-1.) * self.__amplitude
801              for i in xrange(DIM-1):              for i in xrange(DIM-1):
802             x_i=x[i]                 x_i=x[i]
803             min_x=inf(x_i)                 min_x=inf(x_i)
804             max_x=sup(x_i)                 max_x=sup(x_i)
805             k*= sin(self.__n_length*np.pi*(x_i-min_x)/(max_x-min_x))                 k*= sin(self.__n_length*np.pi*(x_i-min_x)/(max_x-min_x))
806              self._reference_data= k              self._reference_data= k
807          return self._reference_data          return self._reference_data
808    
   

Legend:
Removed from v.4125  
changed lines
  Added in v.4132

  ViewVC Help
Powered by ViewVC 1.1.26