feature support specfem3d globe (#143)

* added specfem3d_globe class based on manual inspection of a specfem3d_globe working directory and adjoint simulation

* included laplacian smoothing wrapper and allow user choice between gaussian and laplacian smoothing

* adding specfem3D_globe compatability to SPECFEM Model class. Still need to figure out the intricacies of dealing with multiple regions and whether or not they are updated, which will affect the available kernels, but model reading for 3D globe as well as flavor guessing has been added

* finished adding i/o capabilities for SPECFEM3D_GLOBE models in the Model class. Now allows for user-definition of region for building model, and can read and write globe models which have a different format from specfem2d and 3d models

* bugfix specfem initialize solver ignore directories in DATA directory as we assume we only need text files

* specfem Model class check function can now handle specfem3d_globe models, also deals with all anisotropic model parameters

* specfem3d_globe can also have non-isotropic parameters, making this more general in the check statement

* adjusted specfem solver class (and subclasses) to include a  parameter which is 3dglobe only. this is used to deal with different regions in 3d globe models

* changed specfem3d_globe parent to specfem due to some larger inconsistencies between specfem3d and 3d_globe that make it difficult to have 3d as parent to 3d_globe. need to make some larger scale changes to incorporate

* preprocess read_ascii now deals with the different ascii formats that 3d_globe and 2d/3d have

* bugfix default preprocessing was not properly renaming adjoint sources in line with how specfem3d_globe expects them

* fixed bug for xcombine_sem and xsmooth functions not requiring 'reg' tag for parameter input

seisflows/preprocess/default.py

@@ -329,16 +329,24 @@ def _rename_as_adjoint_source(self, fid):
         Rename synthetic waveforms into filenames consistent with how SPECFEM
         expects adjoint sources to be named. Usually this just means adding
         a '.adj' to the end of the filename
-
-        TODO how does SPECFEM3D_GLOBE expect this? filenames end with .sem.ascii
-            so the .ascii will get replaced. Is that okay?
         """
         if not fid.endswith(".adj"):
             if self.data_format.upper() == "SU":
                 fid = f"{fid}.adj"
             elif self.data_format.upper() == "ASCII":
-                og_extension = os.path.splitext(fid)[-1]  # e.g., .semd
-                fid = fid.replace(og_extension, ".adj")
+                # Differentiate between SPECFEM3D and 3D_GLOBE
+                # SPECFEM3D: NN.SSSS.CCC.sem?
+                # SPECFEM3D_GLOBE: NN.SSSS.CCC.sem.ascii
+                ext = os.path.splitext(fid)[-1]  
+                # SPECFEM3D
+                if ".sem" in ext:
+                    fid = fid.replace(ext, ".adj")
+                # GLOBE (!!! Hardcoded to only work with ASCII format)
+                elif ext == ".ascii":
+                    root, ext1 = os.path.splitext(fid)  # .ascii
+                    root, ext2 = os.path.splitext(root)  # .sem
+                    fid = fid.replace(f"{ext2}{ext1}", ".adj")
+
         return fid
 
     def _setup_quantify_misfit(self, source_name):
@@ -591,12 +599,13 @@ def read_ascii(fid, origintime=None):
     # Honor that Specfem doesn't start exactly on 0
     origintime += times[0]
 
-    # Write out the header information
-    net, sta, cha, fmt = os.path.basename(fid).split('.')
+    # Write out the header information. Deal with the fact that SPECFEM2D/3D and
+    # 3D_GLOBE have slightly different formats for their filenames
+    net, sta, cha, *fmt = os.path.basename(fid).split('.')
     stats = {"network": net, "station": sta, "location": "",
              "channel": cha, "starttime": origintime, "npts": len(data),
              "delta": delta, "mseed": {"dataquality": 'D'},
-             "time_offset": times[0], "format": fmt
+             "time_offset": times[0], "format": fmt[0]
              }
     st = Stream([Trace(data=data, header=stats)])
 

seisflows/solver/specfem.py

@@ -56,10 +56,14 @@ class Specfem:
     :param smooth_h: Gaussian half-width for horizontal smoothing in units
         of meters. If 0., no smoothing applied. Only applicable for workflows:
         ['migration', 'inversion'], ignored for 'forward' workflow.
+        SPECFEM3D_GLOBE only: if `smooth_type`=='laplacian' then this is just 
+        the X and Y extent of the applied smoothing
     :type smooth_h: float
     :param smooth_v: Gaussian half-width for vertical smoothing in units
         of meters. Only applicable for workflows: ['migration', 'inversion'],
         ignored for 'forward' workflow.
+        SPECFEM3D_GLOBE only: if `smooth_type`=='laplacian' then this is just 
+        the Z extent of the applied smoothing
     :type components: str
     :param components: components to consider and tag data with. Should be
         string of letters such as 'RTZ'
@@ -127,7 +131,6 @@ def __init__(self, data_format="ascii",  materials="acoustic",
         self.smooth_h = smooth_h
         self.smooth_v = smooth_v
         self.components = components
-        # self.solver_io = solver_io  # currently not used
         self.source_prefix = source_prefix or "SOURCE"
 
         # Define internally used directory structure
@@ -151,18 +154,22 @@ def __init__(self, data_format="ascii",  materials="acoustic",
 
         self._mpiexec = mpiexec
         self._source_names = None  # for property source_names
-        self._ext = None  # for database file extensions
+        self._ext = ""  # for database file extensions
 
         # Define available choices for check parameters
         self._available_model_types = ["gll"]
         self._available_materials = [
             "ELASTIC", "ACOUSTIC",  # specfem2d, specfem3d
             "ISOTROPIC", "ANISOTROPIC"  # specfem3d_globe
         ]
+        # SPECFEM2D specific attributes. Should be overwritten by 3D versions
         self._available_data_formats = ["ASCII", "SU"]
         self._required_binaries = ["xspecfem2D", "xmeshfem2D", "xcombine_sem",
                                    "xsmooth_sem"]
         self._acceptable_source_prefixes = ["SOURCE", "FORCE", "FORCESOLUTION"]
+
+        # Empty variable that will need to be overwritten by SPECFEM3D_GLOBE
+        self._regions = None
 
     def check(self):
         """
@@ -244,6 +251,21 @@ def check(self):
         assert(isinstance(self.density, bool)), \
             f"solver `density` must be True (variable) or False (constant)"
 
+        # Check the size of the DATA/ directory and let the User know if 
+        # large files are present, e.g., tomo xyz files or topo/bathy
+        for root, dirs, files in os.walk(self.path.specfem_data):
+            for name in files:
+                fullpath = os.path.join(root, name)
+                if not os.path.islink(fullpath):
+                    filesize = os.path.getsize(fullpath) / 1E9  # Bytes -> GB
+                    if filesize > 0.5:
+                        logger.warning(
+                            f"SPECFEM DATA/ file '{fullpath}' is >.5GB and "
+                            f"will be copied {self.ntask} time(s). Please be "
+                            f"sure to check if this file is necessary for your "
+                            f"workflow"
+                            )
+
     @property
     def source_names(self):
         """
@@ -312,6 +334,31 @@ def data_wildcard(self, comp="?"):
         elif self.data_format.upper() == "ASCII":
             return f"*.?X{comp}.sem?"
 
+    def model_wildcard(self, par="*", kernel=False):
+        """
+        Returns a wildcard identifier to search for models kernels generated by
+        the solver. An example SPECFEM2D/3D kernel filename (in 
+        FORTRAN binary file format) is: 'proc000001_rho_kernel.bin'
+        Whereas the corresponding model would be 'proc000001_rho.bin'
+
+        Allows dynamically searching for specific files when renaming, moving
+        or copying files. Also allows for different wildcard for 3D_GLOBE 
+        version
+
+        :type comp: str
+        :param comp: component formatter, defaults to wildcard '?'
+        :type kernel: bool
+        :param kernel: wildcarding a kernel file. If True, adds the 'kernel' 
+            tag. If not, assuming we are wildcarding for a model file
+        :rtype: str
+        :return: wildcard identifier for channels
+        """
+        if kernel:
+            _ker = "_kernel"
+        else:
+            _ker = ""
+        return f"proc??????_{par}{_ker}{self._ext}"
+
     def data_filenames(self, choice="obs"):
         """
         Returns the filenames of SPECFEM2D data, either by the requested
@@ -378,10 +425,11 @@ def model_files(self):
             list of solver parameters
         """
         _model_files = []
-        for parameter in self._parameters:
+        for par in self._parameters:
             _model_files += glob(os.path.join(self.path.mainsolver,
                                               self.model_databases,
-                                              f"*{parameter}{self._ext}"))
+                                              self.model_wildcard(par=par))
+                                              )
         return _model_files
 
     @property
@@ -526,28 +574,32 @@ def adjoint_simulation(self, executables=None, save_kernels=False,
 
         # Rename kernels to work w/ conflicting name conventions
         # Change directory so that the rename doesn't affect the full path
+        # Deals with both SPECFEM3D and 3D_GLOBE, which adds in the 'reg?' tag
         unix.cd(self.kernel_databases)
-        logger.debug(f"renaming output event kernels: 'alpha' -> 'vp'")
-        for tag in ["alpha", "alpha[hv]", "reg1_alpha", "reg1_alpha[hv]"]:
-            names = glob(f"*proc??????_{tag}_kernel{self._ext}")
-            unix.rename(old="alpha", new="vp", names=names)
-
-        logger.debug(f"renaming output event kernels: 'beta' -> 'vs'")
-        for tag in ["beta", "beta[hv]", "reg1_beta", "reg1_beta[hv]"]:
-            names = glob(f"*proc??????_{tag}_kernel{self._ext}")
-            unix.rename(old="beta", new="vs", names=names)
+        for tag in ["alpha", "alpha[hv]", "reg?_alpha", "reg?_alpha[hv]"]:
+            names = glob(self.model_wildcard(par=tag, kernel=True))
+            if names:
+                logger.debug(f"renaming output event kernels: '{tag}' -> 'vp'")
+                unix.rename(old="alpha", new="vp", names=names)
+
+        for tag in ["beta", "beta[hv]", "reg?_beta", "reg?_beta[hv]"]:
+            names = glob(self.model_wildcard(par=tag, kernel=True))
+            if names:
+                logger.debug(f"renaming output event kernels: '{tag}' -> 'vs'")
+                unix.rename(old="beta", new="vs", names=names)
 
         # Save and export the kernels to user-defined locations
         if export_kernels:
             unix.mkdir(export_kernels)
             for par in self._parameters:
-                unix.cp(src=glob(f"*{par}_kernel{self._ext}"),
+                unix.cp(src=glob(self.model_wildcard(par=par, kernel=True)),
                         dst=export_kernels)
 
         if save_kernels:
             unix.mkdir(save_kernels)
             for par in self._parameters:
-                unix.mv(src=glob(f"*{par}_kernel{self._ext}"), dst=save_kernels)
+                unix.mv(src=glob(self.model_wildcard(par=par, kernel=True)),
+                        dst=save_kernels)
 
     def combine(self, input_path, output_path, parameters=None):
         """

seisflows/solver/specfem3d.py

@@ -60,6 +60,7 @@ def __init__(self, source_prefix="CMTSOLUTION", export_vtk=True,
                                    "xgenerate_databases", "xcombine_sem",
                                    "xsmooth_sem", "xcombine_vol_data_vtk"]
 
+        # Internally used parameters set by functions within class
         self._model_databases = None
         self.path._vtk_files = os.path.join(self.path.scratch, "vtk_files")