diff --git a/simpeg_drivers-assets/uijson/fdem1d_forward.ui.json b/simpeg_drivers-assets/uijson/fdem1d_forward.ui.json new file mode 100644 index 00000000..0231c1c2 --- /dev/null +++ b/simpeg_drivers-assets/uijson/fdem1d_forward.ui.json @@ -0,0 +1,212 @@ +{ + "version": "0.3.0-alpha.1", + "title": "Frequency-domain EM (FEM) Forward", + "icon": "surveyairborneem", + "documentation": "https://mirageoscience-simpeg-drivers.readthedocs-hosted.com/en/stable/intro.html", + "conda_environment": "simpeg_drivers", + "run_command": "simpeg_drivers.driver", + "geoh5": "", + "monitoring_directory": "", + "workspace_geoh5": "", + "inversion_type": "fdem 1d", + "physical_property": "conductivity", + "forward_only": true, + "data_object": { + "main": true, + "group": "Survey", + "label": "Object", + "meshType": [ + "{b3a47539-0301-4b27-922e-1dde9d882c60}" + ], + "value": "" + }, + "z_real_channel_bool": { + "group": "Survey", + "main": true, + "label": "Z real component", + "value": true + }, + "z_imag_channel_bool": { + "group": "Survey", + "main": true, + "label": "Z imag component", + "value": true + }, + "v_cell_size": { + "min": 0.0, + "group": "1D Mesh", + "main": true, + "enabled": true, + "label": "Vertical core cell size (m)", + "value": 25.0 + }, + "depth_core": { + "min": 0.0, + "group": "1D Mesh", + "main": true, + "enabled": true, + "label": "Depth of core (m)", + "value": 500.0 + }, + "vertical_padding": { + "min": 0.0, + "group": "1D Mesh", + "main": true, + "dependencyType": "disabled", + "label": "Vertical padding (m)", + "value": 1000.0 + }, + "expansion_factor": { + "main": true, + "group": "1D Mesh", + "label": "Expansion factor", + "value": 1.1 + }, + "mesh": { + "group": "Mesh and models", + "main": true, + "label": "Mesh", + "meshType": "{4ea87376-3ece-438b-bf12-3479733ded46}", + "value": "", + "optional": true, + "enabled": false, + "tooltip": "Select a mesh for the inversion." + }, + "model_type": { + "choiceList": [ + "Conductivity (S/m)", + "Resistivity (Ohm-m)" + ], + "main": true, + "group": "Mesh and models", + "label": "Model units", + "tooltip": "Select the units of the model.", + "value": "Conductivity (S/m)" + }, + "starting_model": { + "association": [ + "Cell", + "Vertex" + ], + "dataType": "Float", + "group": "Mesh and models", + "main": true, + "isValue": true, + "parent": "mesh", + "label": "Value(s)", + "property": "", + "value": 0.001 + }, + "topography_object": { + "main": true, + "group": "Topography", + "label": "Topography", + "meshType": [ + "{202c5db1-a56d-4004-9cad-baafd8899406}", + "{6a057fdc-b355-11e3-95be-fd84a7ffcb88}", + "{f26feba3-aded-494b-b9e9-b2bbcbe298e1}", + "{48f5054a-1c5c-4ca4-9048-80f36dc60a06}", + "{b020a277-90e2-4cd7-84d6-612ee3f25051}" + ], + "value": "", + "optional": true, + "enabled": true, + "tooltip": "Select a topography object to define the active cells for inversion." + }, + "topography": { + "association": [ + "Vertex", + "Cell" + ], + "dataType": "Float", + "group": "Topography", + "main": true, + "optional": true, + "enabled": false, + "label": "Elevation channel", + "tooltip": "Set elevation from channel. If not set the topography will be set from the geometry of the selected 'topography' object.", + "parent": "topography_object", + "dependency": "topography_object", + "dependencyType": "enabled", + "value": "", + "verbose": 2 + }, + "active_model": { + "association": "Cell", + "dataType": [ + "Referenced", + "Boolean", + "Integer" + ], + "group": "Topography", + "main": true, + "enabled": false, + "dependency": "topography_object", + "dependencyType": "disabled", + "label": "Active model", + "tooltip": "Provide the active cell boolean model directly if topography not set.", + "parent": "mesh", + "value": "" + }, + "save_sensitivities": false, + "n_cpu": { + "min": 1, + "group": "Compute", + "optional": true, + "enabled": false, + "label": "Number of CPUs", + "value": 1, + "visible": true + }, + "solver_type": { + "choiceList": [ + "Pardiso", + "Mumps" + ], + "group": "Compute", + "label": "Direct solver", + "tooltip": "Direct solver to use for the forward calculations", + "value": "Pardiso" + }, + "tile_spatial": { + "group": "Compute", + "label": "Number of tiles", + "parent": "data_object", + "isValue": true, + "property": "", + "value": 1, + "min": 1, + "max": 1000, + "verbose": 2, + "tooltip": "Splits the objective function into spatial tiles for distributed computation using the Dask library." + }, + "max_chunk_size": { + "min": 0, + "group": "Compute", + "optional": true, + "enabled": true, + "label": "Maximum chunk size (Mb)", + "value": 128, + "verbose": 3, + "visible": false, + "tooltip": "Limit the chunk size used by Dask for distributed computation." + }, + "out_group": { + "label": "SimPEG group", + "value": "", + "groupType": "{55ed3daf-c192-4d4b-a439-60fa987fe2b8}", + "group": "Drag-and-drop options", + "visible": true, + "optional": true, + "enabled": false, + "tooltip": "Optionally set the SimPEG group to which results will be saved." + }, + "generate_sweep": { + "label": "Generate sweep file", + "group": "Python run preferences", + "main": true, + "value": false, + "tooltip": "Generates a file for sweeping parameters instead of running the application." + }, + "distributed_workers": "" +} diff --git a/simpeg_drivers-assets/uijson/fdem1d_inversion.ui.json b/simpeg_drivers-assets/uijson/fdem1d_inversion.ui.json new file mode 100644 index 00000000..a0d7347b --- /dev/null +++ b/simpeg_drivers-assets/uijson/fdem1d_inversion.ui.json @@ -0,0 +1,578 @@ +{ + "version": "0.3.0-alpha.1", + "title": "Frequency-domain EM 1D (FEM) Inversion", + "icon": "surveyairborneem", + "documentation": "https://mirageoscience-simpeg-drivers.readthedocs-hosted.com/en/stable/intro.html", + "conda_environment": "simpeg_drivers", + "run_command": "simpeg_drivers.driver", + "geoh5": "", + "monitoring_directory": "", + "workspace_geoh5": "", + "inversion_type": "fdem 1d", + "physical_property": "conductivity", + "forward_only": false, + "data_object": { + "main": true, + "group": "Data", + "label": "Object", + "meshType": [ + "{b3a47539-0301-4b27-922e-1dde9d882c60}" + ], + "value": "" + }, + "z_real_channel": { + "association": [ + "Cell", + "Vertex" + ], + "dataType": "Float", + "group": "Data", + "dataGroupType": "Multi-element", + "main": true, + "label": "z-real component (ppm)", + "parent": "data_object", + "optional": true, + "enabled": true, + "value": "" + }, + "z_real_uncertainty": { + "association": [ + "Cell", + "Vertex" + ], + "dataType": "Float", + "group": "Data", + "dataGroupType": "Multi-element", + "main": true, + "label": "Uncertainty (ppm)", + "parent": "data_object", + "dependency": "z_real_channel", + "dependencyType": "enabled", + "value": "" + }, + "z_imag_channel": { + "association": [ + "Cell", + "Vertex" + ], + "dataType": "Float", + "group": "Data", + "dataGroupType": "Multi-element", + "main": true, + "label": "z-imag component (ppm)", + "parent": "data_object", + "optional": true, + "enabled": true, + "value": "" + }, + "z_imag_uncertainty": { + "association": [ + "Cell", + "Vertex" + ], + "dataType": "Float", + "group": "Data", + "dataGroupType": "Multi-element", + "main": true, + "label": "Uncertainty (ppm)", + "parent": "data_object", + "dependency": "z_imag_channel", + "dependencyType": "enabled", + "value": "" + }, + "v_cell_size": { + "min": 0.0, + "group": "1D Mesh", + "main": true, + "enabled": true, + "label": "Vertical core cell size (m)", + "value": 25.0 + }, + "depth_core": { + "min": 0.0, + "group": "1D Mesh", + "main": true, + "enabled": true, + "label": "Depth of core (m)", + "value": 500.0 + }, + "vertical_padding": { + "min": 0.0, + "group": "1D Mesh", + "main": true, + "dependencyType": "disabled", + "label": "Vertical padding (m)", + "value": 1000.0 + }, + "expansion_factor": { + "main": true, + "group": "1D Mesh", + "label": "Expansion factor", + "value": 1.1 + }, + "mesh": { + "group": "Mesh and models", + "main": true, + "label": "Mesh", + "meshType": "{4ea87376-3ece-438b-bf12-3479733ded46}", + "value": "", + "optional": true, + "enabled": false, + "tooltip": "Select a mesh for the inversion." + }, + "model_type": { + "choiceList": [ + "Conductivity (S/m)", + "Resistivity (Ohm-m)" + ], + "main": true, + "group": "Mesh and models", + "label": "Model units", + "tooltip": "Select the units of the model.", + "value": "Conductivity (S/m)" + }, + "starting_model": { + "association": [ + "Cell", + "Vertex" + ], + "dataType": "Float", + "group": "Mesh and models", + "main": true, + "isValue": true, + "parent": "mesh", + "label": "Initial", + "property": "", + "value": 0.001 + }, + "reference_model": { + "association": [ + "Cell", + "Vertex" + ], + "dataType": "Float", + "main": true, + "group": "Mesh and models", + "isValue": true, + "parent": "mesh", + "label": "Reference", + "property": "", + "optional": true, + "enabled": false, + "value": 0.001 + }, + "lower_bound": { + "association": [ + "Cell", + "Vertex" + ], + "main": true, + "dataType": "Float", + "group": "Mesh and models", + "isValue": true, + "parent": "mesh", + "label": "Lower bound", + "property": "", + "optional": true, + "value": 1e-08, + "enabled": false + }, + "upper_bound": { + "association": [ + "Cell", + "Vertex" + ], + "main": true, + "dataType": "Float", + "group": "Mesh and models", + "isValue": true, + "parent": "mesh", + "label": "Upper bound", + "property": "", + "optional": true, + "value": 100.0, + "enabled": false + }, + "topography_object": { + "main": true, + "group": "Topography", + "label": "Topography", + "meshType": [ + "{202c5db1-a56d-4004-9cad-baafd8899406}", + "{6a057fdc-b355-11e3-95be-fd84a7ffcb88}", + "{f26feba3-aded-494b-b9e9-b2bbcbe298e1}", + "{48f5054a-1c5c-4ca4-9048-80f36dc60a06}", + "{b020a277-90e2-4cd7-84d6-612ee3f25051}" + ], + "value": "", + "optional": true, + "enabled": true, + "tooltip": "Select a topography object to define the active cells for inversion." + }, + "topography": { + "association": [ + "Vertex", + "Cell" + ], + "dataType": "Float", + "group": "Topography", + "main": true, + "optional": true, + "enabled": false, + "label": "Elevation channel", + "tooltip": "Set elevation from channel. If not set the topography will be set from the geometry of the selected 'topography' object.", + "parent": "topography_object", + "dependency": "topography_object", + "dependencyType": "enabled", + "value": "", + "verbose": 2 + }, + "active_model": { + "association": "Cell", + "dataType": [ + "Referenced", + "Boolean", + "Integer" + ], + "group": "Topography", + "main": true, + "enabled": false, + "dependency": "topography_object", + "dependencyType": "disabled", + "label": "Active model", + "tooltip": "Provide the active cell boolean model directly if topography not set.", + "parent": "mesh", + "value": "" + }, + "alpha_s": { + "min": 0.0, + "group": "Regularization", + "label": "Reference weight", + "value": 1.0, + "tooltip": "Constant ratio compared to other weights. Larger values result in models that remain close to the reference model", + "dependency": "reference_model", + "dependencyType": "enabled", + "isValue": true, + "parent": "mesh", + "association": "Cell", + "dataType": "Float", + "property": "", + "enabled": true + }, + "length_scale_x": { + "min": 0.0, + "group": "Regularization", + "label": "X-smoothness weight", + "tooltip": "Larger values relative to other smoothness weights will result in x biased smoothness", + "value": 1.0, + "isValue": true, + "parent": "mesh", + "association": "Cell", + "dataType": "Float", + "property": "", + "enabled": true + }, + "length_scale_z": { + "min": 0.0, + "group": "Regularization", + "label": "Z-smoothness weight", + "tooltip": "Larger values relative to other smoothness weights will result in z biased smoothess", + "value": 1.0, + "isValue": true, + "parent": "mesh", + "association": "Cell", + "dataType": "Float", + "property": "", + "enabled": true + }, + "s_norm": { + "association": "Cell", + "dataType": "Float", + "group": "Sparse/blocky model", + "label": "Smallness norm", + "isValue": true, + "parent": "mesh", + "property": "", + "value": 0.0, + "min": 0.0, + "max": 2.0, + "precision": 2, + "lineEdit": true, + "enabled": true, + "dependency": "reference_model", + "dependencyType": "enabled", + "tooltip": "Lp-norm used in the smallness term of the objective function." + }, + "x_norm": { + "association": "Cell", + "dataType": "Float", + "group": "Sparse/blocky model", + "label": "X-smoothness norm", + "isValue": true, + "parent": "mesh", + "property": "", + "value": 2.0, + "min": 0.0, + "max": 2.0, + "precision": 2, + "lineEdit": false, + "enabled": true, + "tooltip": "Lp-norm used in the x-smoothness term of the objective function." + }, + "z_norm": { + "association": "Cell", + "dataType": "Float", + "group": "Sparse/blocky model", + "label": "Z-smoothness norm", + "isValue": true, + "parent": "mesh", + "property": "", + "value": 2.0, + "min": 0.0, + "max": 2.0, + "precision": 2, + "lineEdit": false, + "enabled": true, + "tooltip": "Lp-norm used in the z-smoothness term of the objective function." + }, + "gradient_type": { + "choiceList": [ + "total", + "components" + ], + "group": "Sparse/blocky model", + "label": "Gradient type", + "value": "total", + "verbose": 3, + "tooltip": "Choose whether the IRLS weights for regularization terms are updated using total or partial gradients." + }, + "max_irls_iterations": { + "min": 0, + "group": "Sparse/blocky model", + "label": "Maximum IRLS iterations", + "tooltip": "Incomplete Re-weighted Least Squares iterations for non-L2 problems", + "value": 25, + "enabled": true, + "verbose": 2 + }, + "starting_chi_factor": { + "group": "Sparse/blocky model", + "label": "IRLS start chi factor", + "enabled": true, + "value": 1.0, + "tooltip": "This chi factor will be used to determine the misfit threshold after which IRLS iterations begin.", + "verbose": 3 + }, + "beta_tol": { + "group": "Update IRLS directive", + "label": "Beta tolerance", + "value": 0.5, + "min": 0.0001, + "verbose": 3, + "visible": false + }, + "percentile": { + "group": "Update IRLS directive", + "label": "Percentile", + "value": 95, + "max": 100, + "min": 5, + "verbose": 3, + "visible": false + }, + "chi_factor": { + "min": 0.1, + "max": 20.0, + "precision": 1, + "lineEdit": false, + "group": "Cooling schedule/target", + "label": "Chi factor", + "value": 1.0, + "enabled": true, + "tooltip": "The global target data misfit value." + }, + "auto_scale_misfits": { + "group": "Cooling schedule/target", + "label": "Auto-scale misfits", + "value": false, + "verbose": 3, + "visible": false, + "tooltip": "Whether to auto-scale misfits functions (tile, frequency, joint methods) based on chi-factor." + }, + "initial_beta_ratio": { + "min": 0.0, + "precision": 2, + "group": "Cooling schedule/target", + "optional": true, + "enabled": true, + "label": "Initial beta ratio", + "value": 100.0, + "verbose": 2, + "tooltip": "Estimate the trade-off parameter by scaling the ratio between the largest derivatives in the objective function gradients." + }, + "initial_beta": { + "min": 0.0, + "group": "Cooling schedule/target", + "optional": true, + "enabled": false, + "dependency": "initial_beta_ratio", + "dependencyType": "disabled", + "label": "Initial beta", + "value": 1.0, + "verbose": 2, + "tooltip": "Trade-off parameter between data misfit and regularization." + }, + "cooling_factor": { + "group": "Cooling schedule/target", + "label": "Beta cooling factor", + "tooltip": "Each beta cooling step will be calculated by dividing the current beta by this factor.", + "value": 2.0, + "min": 1.1, + "max": 100, + "precision": 1, + "lineEdit": false, + "verbose": 2 + }, + "cooling_rate": { + "group": "Optimization", + "label": "Iterations per beta", + "value": 2, + "min": 1, + "LineEdit": false, + "max": 10, + "precision": 1, + "verbose": 2, + "enabled": true, + "tooltip": "Set the number of iterations per beta value. Use higher values for more non-linear optimization problems." + }, + "epsilon_cooling_factor": 1.2, + "max_global_iterations": { + "min": 1, + "lineEdit": false, + "group": "Optimization", + "label": "Maximum iterations", + "tooltip": "Number of L2 and IRLS iterations combined", + "value": 50, + "enabled": true + }, + "max_line_search_iterations": { + "group": "Optimization", + "label": "Maximum number of line searches", + "value": 20, + "min": 1, + "enabled": true, + "verbose": 3, + "tooltip": "Perform an Armijo backtracking linesearch for the provided number of iterations." + }, + "max_cg_iterations": { + "min": 0, + "group": "Optimization", + "label": "Maximum CG iterations", + "value": 30, + "enabled": true, + "verbose": 2 + }, + "tol_cg": { + "min": 0, + "group": "Optimization", + "label": "Conjugate gradient tolerance", + "value": 0.0001, + "enabled": true, + "verbose": 3 + }, + "f_min_change": { + "group": "Optimization", + "label": "Minimum change in objective function", + "value": 0.01, + "min": 1e-06, + "verbose": 3, + "enabled": true, + "tooltip": "Minimum decrease in regularization beyond which the IRLS procedure is deemed to have completed." + }, + "sens_wts_threshold": { + "group": "Update sensitivity weights directive", + "tooltip": "Update sensitivity weight threshold", + "label": "Threshold (%)", + "value": 100.0, + "max": 100.0, + "min": 0.0, + "precision": 3, + "enabled": true, + "verbose": 2 + }, + "every_iteration_bool": { + "group": "Update sensitivity weights directive", + "tooltip": "Update weights at every iteration", + "label": "Every iteration", + "value": true, + "verbose": 2, + "enabled": true + }, + "save_sensitivities": { + "group": "Update sensitivity weights directive", + "label": "Save sensitivities", + "tooltip": "Save the summed square row sensitivities to geoh5.", + "value": false + }, + "n_cpu": { + "min": 1, + "group": "Compute", + "optional": true, + "enabled": false, + "label": "Number of CPUs", + "value": 1, + "visible": false + }, + "solver_type": { + "choiceList": [ + "Pardiso", + "Mumps" + ], + "group": "Compute", + "label": "Direct solver", + "tooltip": "Direct solver to use for the forward calculations", + "value": "Pardiso" + }, + "tile_spatial": { + "group": "Compute", + "label": "Number of tiles", + "parent": "data_object", + "isValue": true, + "property": "", + "value": 1, + "min": 1, + "max": 1000, + "verbose": 2, + "visible": false, + "tooltip": "Splits the objective function into spatial tiles for distributed computation using the Dask library." + }, + "store_sensitivities": { + "choiceList": [ + "disk", + "ram" + ], + "group": "Compute", + "label": "Storage device", + "tooltip": "Use disk on a fast local SSD, and RAM elsewhere", + "value": "ram", + "visible": false + }, + "out_group": { + "label": "SimPEG group", + "value": "", + "groupType": "{55ed3daf-c192-4d4b-a439-60fa987fe2b8}", + "group": "Drag-and-drop options", + "visible": true, + "optional": true, + "enabled": false, + "tooltip": "Optionally set the SimPEG group to which results will be saved." + }, + "generate_sweep": { + "label": "Generate sweep file", + "group": "Python run preferences", + "main": true, + "value": false, + "tooltip": "Generates a file for sweeping parameters instead of running the application." + }, + "distributed_workers": "" +} diff --git a/simpeg_drivers-assets/uijson/tdem1d_inversion.ui.json b/simpeg_drivers-assets/uijson/tdem1d_inversion.ui.json index 9a01c302..d173d1a1 100644 --- a/simpeg_drivers-assets/uijson/tdem1d_inversion.ui.json +++ b/simpeg_drivers-assets/uijson/tdem1d_inversion.ui.json @@ -531,7 +531,8 @@ "group": "Compute", "label": "Storage device", "tooltip": "Only RAM storage available for now.", - "value": "ram" + "value": "ram", + "visible": false }, "out_group": { "label": "SimPEG group", diff --git a/simpeg_drivers/__init__.py b/simpeg_drivers/__init__.py index bb0f0078..c29ec416 100644 --- a/simpeg_drivers/__init__.py +++ b/simpeg_drivers/__init__.py @@ -89,13 +89,20 @@ def assets_path() -> Path: "simpeg_drivers.joint.joint_surveys.driver", {"inversion": "JointSurveyDriver"}, ), - "fem": ( + "fdem": ( "simpeg_drivers.electromagnetics.frequency_domain.driver", { - "forward": "FrequenceyDomainElectromagneticsForwardDriver", + "forward": "FDEMForwardDriver", "inversion": "FDEMInversionDriver", }, ), + "fdem 1d": ( + "simpeg_drivers.electromagnetics.frequency_domain_1d.driver", + { + "forward": "FDEM1DForwardDriver", + "inversion": "FDEM1DInversionDriver", + }, + ), "joint cross gradient": ( "simpeg_drivers.joint.joint_cross_gradient.driver", {"inversion": "JointCrossGradientDriver"}, diff --git a/simpeg_drivers/components/data.py b/simpeg_drivers/components/data.py index 735eed3b..af65554c 100644 --- a/simpeg_drivers/components/data.py +++ b/simpeg_drivers/components/data.py @@ -202,7 +202,8 @@ def save_data(self): "magnetotellurics", "tipper", "tdem", - "fem", + "fdem", + "fdem 1d", "tdem 1d", ]: for component, channels in data.items(): @@ -309,7 +310,7 @@ def get_normalizations(self): elif self.params.inversion_type in ["tipper"]: if "imag" in comp: normalizations[chan][comp] = -1 * np.ones(self.mask.sum()) - elif self.params.inversion_type in ["fem"]: + elif "fdem" == self.params.inversion_type: # Assume always ppm data mu0 = 4 * np.pi * 1e-7 offsets = self.params.tx_offsets offsets = { diff --git a/simpeg_drivers/components/factories/directives_factory.py b/simpeg_drivers/components/factories/directives_factory.py index 69214a04..c06f35a7 100644 --- a/simpeg_drivers/components/factories/directives_factory.py +++ b/simpeg_drivers/components/factories/directives_factory.py @@ -224,7 +224,7 @@ def save_iteration_residual_directive(self): if ( self._save_iteration_residual_directive is None and self.factory_type - not in ["tdem", "tdem 1d", "fem", "magnetotellurics", "tipper"] + not in ["tdem", "tdem 1d", "fdem", "fdem 1d", "magnetotellurics", "tipper"] ): self._save_iteration_residual_directive = SaveDataGeoh5Factory( self.params @@ -374,7 +374,8 @@ def assemble_keyword_arguments( "tipper", "tdem", "tdem 1d", - "fem", + "fdem", + "fdem 1d", ]: expmap = maps.ExpMap(inversion_object.mesh) kwargs["transforms"] = [ @@ -464,7 +465,8 @@ def assemble_keyword_arguments( name=None, ): if self.factory_type in [ - "fem", + "fdem", + "fdem 1d", "tdem", "tdem 1d", "magnetotellurics", diff --git a/simpeg_drivers/components/factories/misfit_factory.py b/simpeg_drivers/components/factories/misfit_factory.py index f9d88fdc..52c6c4f8 100644 --- a/simpeg_drivers/components/factories/misfit_factory.py +++ b/simpeg_drivers/components/factories/misfit_factory.py @@ -62,7 +62,7 @@ def assemble_arguments( # pylint: disable=arguments-differ active_cells, ): # Base slice over frequencies - if self.factory_type in ["magnetotellurics", "tipper", "fem"]: + if self.factory_type in ["magnetotellurics", "tipper", "fdem"]: channels = np.unique([list(v) for v in inversion_data.observed.values()]) else: channels = [None] @@ -96,7 +96,7 @@ def assemble_arguments( # pylint: disable=arguments-differ if count == 0: if self.factory_type in [ - "fem", + "fdem", "tdem", "magnetotellurics", "tipper", @@ -209,7 +209,7 @@ def create_nested_simulation( padding_cells=padding_cells, ) inv_type = inversion_data.params.inversion_type - if inv_type in ["fem", "tdem"]: + if inv_type in ["fdem", "tdem"]: compute_em_projections(inversion_data, local_sim) elif ("current" in inv_type or "polarization" in inv_type) and ( "2d" not in inv_type or "pseudo" in inv_type diff --git a/simpeg_drivers/components/factories/receiver_factory.py b/simpeg_drivers/components/factories/receiver_factory.py index c40c60b2..a27e9060 100644 --- a/simpeg_drivers/components/factories/receiver_factory.py +++ b/simpeg_drivers/components/factories/receiver_factory.py @@ -60,9 +60,12 @@ def concrete_object(self): return receivers.Dipole - elif "fem" in self.factory_type: + elif "fdem" in self.factory_type: from simpeg.electromagnetics.frequency_domain import receivers + if "1d" in self.factory_type: + return receivers.PointMagneticFieldSecondary + return receivers.PointMagneticFluxDensitySecondary elif "tdem" in self.factory_type: @@ -129,15 +132,18 @@ def assemble_keyword_arguments( else: kwargs["storeProjections"] = True - if self.factory_type in ["fem", "magnetotellurics", "tipper"]: + if self.factory_type in ["fdem", "fdem 1d", "magnetotellurics", "tipper"]: comp = component.split("_")[0] - kwargs["orientation"] = comp[0] if self.factory_type == "fem" else comp[1:] + kwargs["orientation"] = comp[0] if "fdem" in self.factory_type else comp[1:] kwargs["component"] = component.split("_")[1] if self.factory_type in ["tipper"]: kwargs["orientation"] = kwargs["orientation"][::-1] if "tdem" in self.factory_type: kwargs["orientation"] = component + if self.factory_type == "fdem 1d": + kwargs["data_type"] = "ppm" + return kwargs def build(self, locations=None, data=None, local_index=None, component=None): diff --git a/simpeg_drivers/components/factories/simpeg_factory.py b/simpeg_drivers/components/factories/simpeg_factory.py index b9611bdf..44d92db0 100644 --- a/simpeg_drivers/components/factories/simpeg_factory.py +++ b/simpeg_drivers/components/factories/simpeg_factory.py @@ -56,7 +56,8 @@ class SimPEGFactory(ABC): "induced polarization 3d", "induced polarization 2d", "induced polarization pseudo 3d", - "fem", + "fdem", + "fdem 1d", "tdem", "tdem 1d", "magnetotellurics", diff --git a/simpeg_drivers/components/factories/simulation_factory.py b/simpeg_drivers/components/factories/simulation_factory.py index 1da9656b..e2567b09 100644 --- a/simpeg_drivers/components/factories/simulation_factory.py +++ b/simpeg_drivers/components/factories/simulation_factory.py @@ -49,7 +49,7 @@ def __init__(self, params: BaseParams | BaseOptions): "induced polarization pseudo 3d", "magnetotellurics", "tipper", - "fem", + "fdem", "tdem", ]: import pymatsolver.direct as solver_module @@ -97,11 +97,16 @@ def concrete_object(self): return simulation.Simulation3DPrimarySecondary - if self.factory_type in ["fem"]: + if self.factory_type in ["fdem"]: from simpeg.electromagnetics.frequency_domain import simulation return simulation.Simulation3DMagneticFluxDensity + if self.factory_type in ["fdem 1d"]: + from simpeg.electromagnetics.frequency_domain import simulation_1d + + return simulation_1d.Simulation1DLayered + if self.factory_type in ["tdem"]: from simpeg.electromagnetics.time_domain import simulation @@ -179,7 +184,7 @@ def assemble_keyword_arguments( "direct current 2d", "magnetotellurics", "tipper", - "fem", + "fdem", "tdem", ]: actmap = maps.InjectActiveCells( @@ -194,7 +199,7 @@ def assemble_keyword_arguments( * self.params.unit_conversion ) - if self.factory_type in ["tdem 1d"]: + if "1d" in self.factory_type: kwargs["sigmaMap"] = maps.ExpMap(mesh) kwargs["thicknesses"] = local_mesh.h[0][1:][::-1] kwargs["topo"] = active_cells[tile_id] diff --git a/simpeg_drivers/components/factories/source_factory.py b/simpeg_drivers/components/factories/source_factory.py index c52c0e41..15ee7cc5 100644 --- a/simpeg_drivers/components/factories/source_factory.py +++ b/simpeg_drivers/components/factories/source_factory.py @@ -56,23 +56,29 @@ def concrete_object(self): elif "induced polarization" in self.factory_type: return dc_sources.Dipole - elif "fem" in self.factory_type: + elif "fdem" in self.factory_type: + if "fdem 1d" == self.factory_type and np.allclose( + np.kron( + np.ones((len(self.params.data_object.channels), 1)), + self.params.data_object.vertices, + ), + self.params.data_object.complement.vertices, + ): + return fem_sources.CircularLoop + return fem_sources.MagDipole - elif "tdem" == self.factory_type: + elif "tdem" in self.factory_type: if isinstance(self.params.data_object, LargeLoopGroundTEMReceivers): return tem_sources.LineCurrent - else: - return tem_sources.MagDipole - elif "tdem 1d" == self.factory_type: - if np.allclose( + if "tdem 1d" == self.factory_type and np.allclose( self.params.data_object.vertices, self.params.data_object.complement.vertices, ): return tem_sources.CircularLoop - else: - return tem_sources.MagDipole + + return tem_sources.MagDipole elif self.factory_type in ["magnetotellurics", "tipper"]: return ns_sources.PlanewaveXYPrimary @@ -101,7 +107,7 @@ def assemble_arguments( locations=locations, ) - elif self.factory_type in ["fem", "magnetotellurics", "tipper"]: + elif self.factory_type in ["fdem", "fdem 1d", "magnetotellurics", "tipper"]: args.append(receivers) args.append(frequency) @@ -133,16 +139,16 @@ def assemble_keyword_arguments( # pylint: disable=arguments-differ kwargs["sigma_primary"] = [background] - if self.factory_type in ["fem"]: + if "fdem" in self.factory_type: kwargs["location"] = locations if "tdem" in self.factory_type: kwargs["location"] = locations kwargs["waveform"] = waveform - if self.factory_type == "tdem 1d": + if "1d" in self.factory_type: if isinstance( self.concrete_object(), - tem_sources.CircularLoop, + tem_sources.CircularLoop | fem_sources.CircularLoop, ): kwargs["moment"] = 1.0 else: diff --git a/simpeg_drivers/components/factories/survey_factory.py b/simpeg_drivers/components/factories/survey_factory.py index 9abdc1ad..9c9a5fb1 100644 --- a/simpeg_drivers/components/factories/survey_factory.py +++ b/simpeg_drivers/components/factories/survey_factory.py @@ -94,7 +94,7 @@ def concrete_object(self): elif "induced polarization" in self.factory_type: from simpeg.electromagnetics.static.induced_polarization import survey - elif "fem" in self.factory_type: + elif "fdem" in self.factory_type: from simpeg.electromagnetics.frequency_domain import survey elif "tdem" in self.factory_type: @@ -124,11 +124,11 @@ def assemble_arguments(self, data=None, local_index=None, channel=None): if "current" in self.factory_type or "polarization" in self.factory_type: return self._dcip_arguments(data=data, local_index=local_index) - elif self.factory_type in ["tdem", "tdem 1d"]: + elif "tdem" in self.factory_type: return self._tdem_arguments(data=data) elif self.factory_type in ["magnetotellurics", "tipper"]: return self._naturalsource_arguments(data=data, frequency=channel) - elif self.factory_type in ["fem"]: + elif "fdem" in self.factory_type: return self._fem_arguments(data=data, channel=channel) else: receivers = ReceiversFactory(self.params).build( @@ -194,7 +194,7 @@ def _add_data(self, survey, data, local_index, channel): if isinstance(local_index, list): local_index = np.hstack(local_index) - if self.factory_type in ["fem", "tdem", "tdem 1d"]: + if self.factory_type in ["fdem", "fdem 1d", "tdem", "tdem 1d"]: dobs = [] uncerts = [] diff --git a/simpeg_drivers/driver.py b/simpeg_drivers/driver.py index c80a8ee4..4fe17c2e 100644 --- a/simpeg_drivers/driver.py +++ b/simpeg_drivers/driver.py @@ -449,8 +449,12 @@ def get_regularization(self): ) # Adjustment for 2D versus 3D problems - comps = "sxz" if "2d" in self.params.inversion_type else "sxyz" - avg_comps = "sxy" if "2d" in self.params.inversion_type else "sxyz" + # TODO check this part + is_2d_reg = ( + "2d" in self.params.inversion_type or "1d" in self.params.inversion_type + ) + comps = "sxz" if is_2d_reg else "sxyz" + avg_comps = "sxy" if is_2d_reg else "sxyz" weights = ["alpha_s"] + [f"length_scale_{k}" for k in comps[1:]] for comp, avg_comp, objfct, weight in zip( comps, avg_comps, reg.objfcts, weights diff --git a/simpeg_drivers/electromagnetics/base_1d_driver.py b/simpeg_drivers/electromagnetics/base_1d_driver.py new file mode 100644 index 00000000..96242169 --- /dev/null +++ b/simpeg_drivers/electromagnetics/base_1d_driver.py @@ -0,0 +1,119 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + + +from __future__ import annotations + +from logging import getLogger + +import numpy as np +from discretize import TensorMesh +from discretize.utils import mesh_utils +from geoh5py import Workspace +from geoh5py.shared.merging.drape_model import DrapeModelMerger +from geoh5py.ui_json.ui_json import fetch_active_workspace + +from simpeg_drivers.components.factories import MisfitFactory +from simpeg_drivers.components.meshes import InversionMesh +from simpeg_drivers.driver import InversionDriver +from simpeg_drivers.utils.utils import topo_drape_elevation, xyz_2_drape_model + + +logger = getLogger(__name__) + + +class Base1DDriver(InversionDriver): + """Base 1D driver for electromagnetic simulations.""" + + _params_class = None + _validations = None + + def __init__(self, workspace: Workspace, **kwargs): + super().__init__(workspace, **kwargs) + + self.layers_mesh: TensorMesh = self.get_1d_mesh() + self.topo_z_drape = topo_drape_elevation( + self.params.data_object.vertices, + self.inversion_topography.locations, + ) + + @property + def inversion_mesh(self) -> InversionMesh: + """Inversion mesh""" + if getattr(self, "_inversion_mesh", None) is None: + with fetch_active_workspace(self.workspace, mode="r+"): + drape_models = [] + temp_work = Workspace() + for part in self.params.data_object.unique_parts: + indices = self.params.data_object.parts == part + drape_models.append( + xyz_2_drape_model( + temp_work, + self.topo_z_drape[indices], + self.layers_mesh.h[0][::-1], + ) + ) + + entity = DrapeModelMerger.create_object( + self.workspace, drape_models, parent=self.out_group + ) + + self._inversion_mesh = InversionMesh( + self.workspace, self.params, entity=entity + ) + self._inversion_mesh.layers_mesh = self.layers_mesh + + return self._inversion_mesh + + def get_1d_mesh(self) -> TensorMesh: + layers_mesh = mesh_utils.mesh_builder_xyz( + np.c_[0], + np.r_[self.params.drape_model.v_cell_size], + padding_distance=[ + [self.params.drape_model.vertical_padding, 0], + ], + depth_core=self.params.drape_model.depth_core, + expansion_factor=self.params.drape_model.expansion_factor, + mesh_type="tensor", + ) + return layers_mesh + + @property + def data_misfit(self): + """The Simpeg.data_misfit class""" + if getattr(self, "_data_misfit", None) is None: + with fetch_active_workspace(self.workspace, mode="r+"): + # Tile locations + tiles = self.get_tiles() + + logger.info("Setting up %i tile(s) . . .", len(tiles)) + # Build tiled misfits and combine to form global misfit + self._data_misfit, self._sorting, self._ordering = MisfitFactory( + self.params, models=self.models + ).build( + tiles, + self.inversion_data, + self.inversion_mesh, + self.topo_z_drape, + ) + self.models.active_cells = np.ones( + self.inversion_mesh.mesh.n_cells, dtype=bool + ) + logger.info("Done.") + + self.inversion_data.save_data() + self._data_misfit.multipliers = np.asarray( + self._data_misfit.multipliers, dtype=float + ) + + if self.client: + self.distributed_misfits() + + return self._data_misfit diff --git a/simpeg_drivers/electromagnetics/frequency_domain/params.py b/simpeg_drivers/electromagnetics/frequency_domain/params.py index 1012f6da..f0ae20bf 100644 --- a/simpeg_drivers/electromagnetics/frequency_domain/params.py +++ b/simpeg_drivers/electromagnetics/frequency_domain/params.py @@ -30,25 +30,13 @@ ) -class FDEMForwardOptions(EMDataMixin, BaseForwardOptions): +class BaseFDEMOptions(EMDataMixin): """ - Frequency Domain Electromagnetic Forward options. - - :param z_real_channel_bool: Real impedance channel boolean. - :param z_imag_channel_bool: Imaginary impedance channel boolean. - :param model_type: Specify whether the models are provided in resistivity or conductivity. + Base Frequency Domain Electromagnetic options. """ - name: ClassVar[str] = "Frequency Domain Electromagnetics Forward" - default_ui_json: ClassVar[Path] = assets_path() / "uijson/fem_forward.ui.json" - - title: str = "Frequency-domain EM (FEM) Forward" physical_property: str = "conductivity" - inversion_type: str = "fem" - data_object: Receivers - z_real_channel_bool: bool - z_imag_channel_bool: bool model_type: str = "Conductivity (S/m)" @property @@ -71,12 +59,32 @@ def tx_offsets(self): def unit_conversion(self): """Return time unit conversion factor.""" conversion = { - "Hertz (Hz)": 1.0, + "Seconds (s)": 1.0, + "Milliseconds (ms)": 1e-3, + "Microseconds (us)": 1e-6, } return conversion[self.data_object.unit] -class FDEMInversionOptions(EMDataMixin, BaseInversionOptions): +class FDEMForwardOptions(BaseFDEMOptions, BaseForwardOptions): + """ + Frequency Domain Electromagnetic Forward options. + + :param z_real_channel_bool: Real impedance channel boolean. + :param z_imag_channel_bool: Imaginary impedance channel boolean. + :param model_type: Specify whether the models are provided in resistivity or conductivity. + """ + + name: ClassVar[str] = "Frequency Domain Electromagnetics Forward" + default_ui_json: ClassVar[Path] = assets_path() / "uijson/fem_forward.ui.json" + inversion_type: str = "fdem" + title: str = "Frequency-domain EM (FEM) Forward" + + z_real_channel_bool: bool + z_imag_channel_bool: bool + + +class FDEMInversionOptions(BaseFDEMOptions, BaseInversionOptions): """ Frequency Domain Electromagnetic Inversion options. @@ -89,38 +97,10 @@ class FDEMInversionOptions(EMDataMixin, BaseInversionOptions): name: ClassVar[str] = "Frequency Domain Electromagnetics Inversion" default_ui_json: ClassVar[Path] = assets_path() / "uijson/fem_inversion.ui.json" - + inversion_type: str = "fdem" title: str = "Frequency-domain EM (FEM) Inversion" - physical_property: str = "conductivity" - inversion_type: str = "fem" - data_object: Receivers z_real_channel: PropertyGroup | None = None z_real_uncertainty: PropertyGroup | None = None z_imag_channel: PropertyGroup | None = None z_imag_uncertainty: PropertyGroup | None = None - model_type: str = "Conductivity (S/m)" - - @property - def tx_offsets(self): - """Return transmitter offsets from frequency metadata""" - - try: - offset_data = self.data_object.metadata["EM Dataset"][ - "Frequency configurations" - ] - tx_offsets = {k["Frequency"]: k["Offset"] for k in offset_data} - - except KeyError as exception: - msg = "Metadata must contain 'Frequency configurations' dictionary with 'Offset' data." - raise KeyError(msg) from exception - - return tx_offsets - - @property - def unit_conversion(self): - """Return time unit conversion factor.""" - conversion = { - "Hertz (Hz)": 1.0, - } - return conversion[self.data_object.unit] diff --git a/simpeg_drivers/electromagnetics/frequency_domain_1d/__init__.py b/simpeg_drivers/electromagnetics/frequency_domain_1d/__init__.py new file mode 100644 index 00000000..4d06f672 --- /dev/null +++ b/simpeg_drivers/electromagnetics/frequency_domain_1d/__init__.py @@ -0,0 +1,9 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' diff --git a/simpeg_drivers/electromagnetics/frequency_domain_1d/driver.py b/simpeg_drivers/electromagnetics/frequency_domain_1d/driver.py new file mode 100644 index 00000000..f0e99923 --- /dev/null +++ b/simpeg_drivers/electromagnetics/frequency_domain_1d/driver.py @@ -0,0 +1,33 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + + +from __future__ import annotations + +from simpeg_drivers.electromagnetics.base_1d_driver import Base1DDriver + +from .params import ( + FDEM1DForwardOptions, + FDEM1DInversionOptions, +) + + +class FDEM1DForwardDriver(Base1DDriver): + """Frequency Domain 1D Electromagnetic forward driver.""" + + _params_class = FDEM1DForwardOptions + _validations = None + + +class FDEM1DInversionDriver(Base1DDriver): + """Frequency Domain 1D Electromagnetic inversion driver.""" + + _params_class = FDEM1DInversionOptions + _validations = None diff --git a/simpeg_drivers/electromagnetics/frequency_domain_1d/params.py b/simpeg_drivers/electromagnetics/frequency_domain_1d/params.py new file mode 100644 index 00000000..d56bc2a8 --- /dev/null +++ b/simpeg_drivers/electromagnetics/frequency_domain_1d/params.py @@ -0,0 +1,80 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + + +from __future__ import annotations + +from pathlib import Path +from typing import ClassVar + +from geoh5py.groups import PropertyGroup + +from simpeg_drivers import assets_path +from simpeg_drivers.electromagnetics.frequency_domain import ( + FDEMForwardOptions, + FDEMInversionOptions, +) +from simpeg_drivers.params import DrapeModelOptions + + +class FDEM1DForwardOptions(FDEMForwardOptions): + """ + Frequency Domain Electromagnetic forward options. + + :param z_channel_bool: Z-component data channel boolean. + :param x_channel_bool: X-component data channel boolean. + :param y_channel_bool: Y-component data channel boolean. + :param model_type: Specify whether the models are provided in resistivity or conductivity. + :param data_units: Units for the FEM data + """ + + name: ClassVar[str] = "Frequency Domain 1D Electromagnetics Forward" + default_ui_json: ClassVar[Path] = assets_path() / "uijson/fdem1d_forward.ui.json" + + title: str = "Frequency-domain EM-1D (FEM-1D) Forward" + inversion_type: str = "fdem 1d" + drape_model: DrapeModelOptions = DrapeModelOptions( + u_cell_size=10.0, + v_cell_size=10.0, + depth_core=100.0, + horizontal_padding=0.0, + vertical_padding=100.0, + expansion_factor=1.1, + ) + + +class FDEM1DInversionOptions(FDEMInversionOptions): + """ + Frequency Domain Electromagnetic Inversion options. + + :param z_real_channel: Real Z-component data channel. + :param z_real_uncertainty: Real Z-component data channel uncertainty. + :param z_imag_channel: Imaginary Z-component data channel. + :param z_imag_uncertainty: Imaginary Z-component data channel uncertainty. + :param model_type: Specify whether the models are provided in resistivity or conductivity. + :param data_units: Units for the FEM data + """ + + name: ClassVar[str] = "Frequency Domain 1D Electromagnetics Inversion" + default_ui_json: ClassVar[Path] = assets_path() / "uijson/fdem1d_inversion.ui.json" + + title: str = "Frequency-domain EM-1D (FEM-1D) Inversion" + inversion_type: str = "fdem 1d" + y_norm: None = None + drape_model: DrapeModelOptions = DrapeModelOptions( + u_cell_size=10.0, + v_cell_size=10.0, + depth_core=100.0, + horizontal_padding=0.0, + vertical_padding=100.0, + expansion_factor=1.1, + ) + auto_scale_misfits: bool = False + sens_wts_threshold: float = 100.0 diff --git a/simpeg_drivers/electromagnetics/time_domain/__init__.py b/simpeg_drivers/electromagnetics/time_domain/__init__.py index 54036b7b..4d06f672 100644 --- a/simpeg_drivers/electromagnetics/time_domain/__init__.py +++ b/simpeg_drivers/electromagnetics/time_domain/__init__.py @@ -7,12 +7,3 @@ # (see LICENSE file at the root of this source code package). ' # ' # ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' - - -from .params import ( - TDEMForwardOptions, - TDEMInversionOptions, -) - -# pylint: disable=unused-import -# flake8: noqa diff --git a/simpeg_drivers/electromagnetics/time_domain_1d/driver.py b/simpeg_drivers/electromagnetics/time_domain_1d/driver.py index 1daac739..7727f835 100644 --- a/simpeg_drivers/electromagnetics/time_domain_1d/driver.py +++ b/simpeg_drivers/electromagnetics/time_domain_1d/driver.py @@ -11,17 +11,7 @@ from __future__ import annotations -import numpy as np -from discretize import TensorMesh -from discretize.utils import mesh_utils -from geoh5py import Workspace -from geoh5py.shared.merging.drape_model import DrapeModelMerger -from geoh5py.ui_json.ui_json import fetch_active_workspace - -from simpeg_drivers.components.factories import MisfitFactory -from simpeg_drivers.components.meshes import InversionMesh -from simpeg_drivers.driver import InversionDriver -from simpeg_drivers.utils.utils import topo_drape_elevation, xyz_2_drape_model +from simpeg_drivers.electromagnetics.base_1d_driver import Base1DDriver from .params import ( TDEM1DForwardOptions, @@ -29,106 +19,15 @@ ) -class Base1DDriver(InversionDriver): - """Base 1D driver for electromagnetic simulations.""" - - _params_class = None - _validations = None - - def __init__(self, workspace: Workspace, **kwargs): - super().__init__(workspace, **kwargs) - - self.layers_mesh: TensorMesh = self.get_1d_mesh() - self.topo_z_drape = topo_drape_elevation( - self.params.data_object.vertices, - self.inversion_topography.locations, - ) - - @property - def inversion_mesh(self) -> InversionMesh: - """Inversion mesh""" - if getattr(self, "_inversion_mesh", None) is None: - # temp_workspace = Workspace() - with fetch_active_workspace(self.workspace, mode="r+"): - drape_models = [] - temp_work = Workspace() - for part in self.params.data_object.unique_parts: - indices = self.params.data_object.parts == part - drape_models.append( - xyz_2_drape_model( - temp_work, - self.topo_z_drape[indices], - self.layers_mesh.h[0][::-1], - ) - ) - - entity = DrapeModelMerger.create_object( - self.workspace, drape_models, parent=self.out_group - ) - - self._inversion_mesh = InversionMesh( - self.workspace, self.params, entity=entity - ) - self._inversion_mesh.layers_mesh = self.layers_mesh - - return self._inversion_mesh - - def get_1d_mesh(self) -> TensorMesh: - layers_mesh = mesh_utils.mesh_builder_xyz( - np.c_[0], - np.r_[self.params.drape_model.v_cell_size], - padding_distance=[ - [self.params.drape_model.vertical_padding, 0], - ], - depth_core=self.params.drape_model.depth_core, - expansion_factor=self.params.drape_model.expansion_factor, - mesh_type="tensor", - ) - return layers_mesh - - @property - def data_misfit(self): - """The Simpeg.data_misfit class""" - if getattr(self, "_data_misfit", None) is None: - with fetch_active_workspace(self.workspace, mode="r+"): - # Tile locations - tiles = self.get_tiles() - - print(f"Setting up {len(tiles)} tile(s) . . .") - # Build tiled misfits and combine to form global misfit - self._data_misfit, self._sorting, self._ordering = MisfitFactory( - self.params, models=self.models - ).build( - tiles, - self.inversion_data, - self.inversion_mesh, - self.topo_z_drape, - ) - self.models.active_cells = np.ones( - self.inversion_mesh.mesh.n_cells, dtype=bool - ) - print("Done.") - - self.inversion_data.save_data() - self._data_misfit.multipliers = np.asarray( - self._data_misfit.multipliers, dtype=float - ) - - if self.client: - self.distributed_misfits() - - return self._data_misfit - - class TDEM1DForwardDriver(Base1DDriver): - """Time Domain Electromagnetic forward driver.""" + """Time Domain 1D Electromagnetic forward driver.""" _params_class = TDEM1DForwardOptions _validations = None class TDEM1DInversionDriver(Base1DDriver): - """Time Domain Electromagnetic inversion driver.""" + """Time Domain 1D Electromagnetic inversion driver.""" _params_class = TDEM1DInversionOptions _validations = None diff --git a/simpeg_drivers/electromagnetics/time_domain_1d/params.py b/simpeg_drivers/electromagnetics/time_domain_1d/params.py index f372cb78..cd2bd76a 100644 --- a/simpeg_drivers/electromagnetics/time_domain_1d/params.py +++ b/simpeg_drivers/electromagnetics/time_domain_1d/params.py @@ -54,7 +54,7 @@ class TDEM1DForwardOptions(EMDataMixin, BaseForwardOptions): physical_property: str = "conductivity" data_object: Receivers - z_channel_bool: bool | None = None + z_channel_bool: bool x_channel_bool: None = None y_channel_bool: None = None data_units: str = "dB/dt (T/s)" diff --git a/simpeg_drivers/params.py b/simpeg_drivers/params.py index 41747ba8..9cb41408 100644 --- a/simpeg_drivers/params.py +++ b/simpeg_drivers/params.py @@ -237,7 +237,7 @@ def padding_cells(self) -> int: if self.tile_spatial == 1: return 100 - return 4 if self.inversion_type in ["fem", "tdem"] else 6 + return 4 if self.inversion_type in ["fdem", "tdem"] else 6 class BaseForwardOptions(CoreOptions): diff --git a/simpeg_drivers/utils/testing.py b/simpeg_drivers/utils/testing.py index a171be44..9cc7589d 100644 --- a/simpeg_drivers/utils/testing.py +++ b/simpeg_drivers/utils/testing.py @@ -280,7 +280,7 @@ def topo_drape(x, y): # survey.cells = survey.cells[dist < 100.0, :] survey.remove_cells(np.where(dist > 200)[0]) - elif inversion_type == "fem": + elif "fdem" in inversion_type: survey = AirborneFEMReceivers.create( geoh5, vertices=vertices, name="Airborne_rx" ) @@ -466,7 +466,7 @@ def topo_drape(x, y): finalize=False, ) - if inversion_type in ["fem", "airborne_tem"]: + if inversion_type in ["fdem", "airborne_tem"]: mesh = OctreeDriver.refine_tree_from_points( mesh, vertices, diff --git a/tests/run_tests/driver_airborne_fem_1d_test.py b/tests/run_tests/driver_airborne_fem_1d_test.py new file mode 100644 index 00000000..12e090be --- /dev/null +++ b/tests/run_tests/driver_airborne_fem_1d_test.py @@ -0,0 +1,182 @@ +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' +# Copyright (c) 2025 Mira Geoscience Ltd. ' +# ' +# This file is part of simpeg-drivers package. ' +# ' +# simpeg-drivers is distributed under the terms and conditions of the MIT License ' +# (see LICENSE file at the root of this source code package). ' +# ' +# ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' + +# pylint: disable=too-many-locals + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +from geoh5py import Workspace +from geoh5py.groups import SimPEGGroup + +from simpeg_drivers.electromagnetics.frequency_domain_1d.driver import ( + FDEM1DForwardDriver, + FDEM1DInversionDriver, +) +from simpeg_drivers.electromagnetics.frequency_domain_1d.params import ( + FDEM1DForwardOptions, + FDEM1DInversionOptions, +) +from simpeg_drivers.params import ActiveCellsOptions +from simpeg_drivers.utils.testing import check_target, setup_inversion_workspace +from simpeg_drivers.utils.utils import get_inversion_output + + +# To test the full run and validate the inversion. +# Move this file out of the test directory and run. + +target_run = {"data_norm": 638.15180, "phi_d": 59500, "phi_m": 168} + + +def test_fem_fwr_1d_run( + tmp_path: Path, + n_grid_points=3, + refinement=(2,), + cell_size=(20.0, 20.0, 20.0), +): + # Run the forward + geoh5, _, model, survey, topography = setup_inversion_workspace( + tmp_path, + background=1e-4, + anomaly=0.1, + n_electrodes=n_grid_points, + n_lines=n_grid_points, + refinement=refinement, + drape_height=10.0, + cell_size=cell_size, + padding_distance=400, + inversion_type="fdem 1d", + flatten=False, + ) + params = FDEM1DForwardOptions( + geoh5=geoh5, + mesh=model.parent, + active_cells=ActiveCellsOptions(topography_object=topography), + data_object=survey, + starting_model=model, + z_real_channel_bool=True, + z_imag_channel_bool=True, + ) + + fwr_driver = FDEM1DForwardDriver(params) + fwr_driver.run() + + +def test_fem_1d_run(tmp_path: Path, max_iterations=1, pytest=True): + workpath = tmp_path / "inversion_test.ui.geoh5" + if pytest: + workpath = tmp_path.parent / "test_fem_fwr_1d_run0" / "inversion_test.ui.geoh5" + + with Workspace(workpath) as geoh5: + survey = next( + child + for child in geoh5.get_entity("Airborne_rx") + if not isinstance(child.parent, SimPEGGroup) + ) + mesh = geoh5.get_entity("mesh")[0] + topography = geoh5.get_entity("topography")[0] + data = {} + uncertainties = {} + components = { + "z_real": "z_real", + "z_imag": "z_imag", + } + + for comp, cname in components.items(): + data[cname] = [] + uncertainties[f"{cname} uncertainties"] = [] + for ind, freq in enumerate(survey.channels): + data_entity = geoh5.get_entity(f"Iteration_0_{comp}_[{ind}]")[0].copy( + parent=survey + ) + data[cname].append(data_entity) + abs_val = np.abs(data_entity.values) + uncert = survey.add_data( + { + f"uncertainty_{comp}_[{ind}]": { + "values": np.ones_like(abs_val) + * freq + / 200.0 # * 2**(np.abs(ind-1)) + } + } + ) + uncertainties[f"{cname} uncertainties"].append( + uncert.copy(parent=survey) + ) + + data_groups = survey.add_components_data(data) + uncert_groups = survey.add_components_data(uncertainties) + + data_kwargs = {} + for comp, data_group, uncert_group in zip( + components, data_groups, uncert_groups, strict=True + ): + data_kwargs[f"{comp}_channel"] = data_group + data_kwargs[f"{comp}_uncertainty"] = uncert_group + + orig_z_real_1 = geoh5.get_entity("Iteration_0_z_real_[0]")[0].values + + # Run the inverse + params = FDEM1DInversionOptions( + geoh5=geoh5, + mesh=mesh, + active_cells=ActiveCellsOptions(topography_object=topography), + data_object=survey, + starting_model=1e-3, + reference_model=1e-3, + alpha_s=0.0, + s_norm=0.0, + x_norm=0.0, + z_norm=0.0, + gradient_type="components", + upper_bound=0.75, + max_global_iterations=max_iterations, + initial_beta_ratio=1e1, + percentile=100, + cooling_rate=3, + chi_factor=0.25, + store_sensitivities="ram", + sens_wts_threshold=1.0, + **data_kwargs, + ) + params.write_ui_json(path=tmp_path / "Inv_run.ui.json") + driver = FDEM1DInversionDriver(params) + driver.run() + + with geoh5.open() as run_ws: + output = get_inversion_output( + driver.params.geoh5.h5file, driver.params.out_group.uid + ) + output["data"] = orig_z_real_1 + + assert ( + run_ws.get_entity("Iteration_1_z_imag_[1]")[0].entity_type.uid + == run_ws.get_entity("Observed_z_imag_[1]")[0].entity_type.uid + ) + + if pytest: + check_target(output, target_run) + nan_ind = np.isnan(run_ws.get_entity("Iteration_0_model")[0].values) + inactive_ind = run_ws.get_entity("active_cells")[0].values == 0 + assert np.all(nan_ind == inactive_ind) + + +if __name__ == "__main__": + # Full run + test_fem_fwr_1d_run( + Path("./"), n_grid_points=5, cell_size=(5.0, 5.0, 5.0), refinement=(4, 4, 4) + ) + test_fem_1d_run( + Path("./"), + max_iterations=15, + pytest=False, + ) diff --git a/tests/run_tests/driver_airborne_tem_1d_test.py b/tests/run_tests/driver_airborne_tem_1d_test.py index e38b6714..452f9044 100644 --- a/tests/run_tests/driver_airborne_tem_1d_test.py +++ b/tests/run_tests/driver_airborne_tem_1d_test.py @@ -32,7 +32,7 @@ # To test the full run and validate the inversion. # Move this file out of the test directory and run. -target_run = {"data_norm": 6.15712e-10, "phi_d": 110, "phi_m": 89800} +target_run = {"data_norm": 6.15712e-10, "phi_d": 109, "phi_m": 102000} def test_airborne_tem_1d_fwr_run( diff --git a/tests/run_tests/driver_airborne_tem_test.py b/tests/run_tests/driver_airborne_tem_test.py index f2b3abc5..b90de30f 100644 --- a/tests/run_tests/driver_airborne_tem_test.py +++ b/tests/run_tests/driver_airborne_tem_test.py @@ -17,14 +17,14 @@ from geoh5py.workspace import Workspace from pytest import raises -from simpeg_drivers.electromagnetics.time_domain import ( - TDEMForwardOptions, - TDEMInversionOptions, -) from simpeg_drivers.electromagnetics.time_domain.driver import ( TDEMForwardDriver, TDEMInversionDriver, ) +from simpeg_drivers.electromagnetics.time_domain.params import ( + TDEMForwardOptions, + TDEMInversionOptions, +) from simpeg_drivers.params import ActiveCellsOptions from simpeg_drivers.utils.testing import check_target, setup_inversion_workspace from simpeg_drivers.utils.utils import get_inversion_output diff --git a/tests/run_tests/driver_fem_test.py b/tests/run_tests/driver_fem_test.py index 1d77f573..b5d95e3c 100644 --- a/tests/run_tests/driver_fem_test.py +++ b/tests/run_tests/driver_fem_test.py @@ -54,7 +54,7 @@ def test_fem_fwr_run( drape_height=15.0, cell_size=cell_size, padding_distance=400, - inversion_type="fem", + inversion_type="fdem", flatten=True, ) params = FDEMForwardOptions( diff --git a/tests/run_tests/driver_ground_tem_test.py b/tests/run_tests/driver_ground_tem_test.py index b5f7fa3c..46a257d4 100644 --- a/tests/run_tests/driver_ground_tem_test.py +++ b/tests/run_tests/driver_ground_tem_test.py @@ -17,14 +17,14 @@ from geoh5py.workspace import Workspace from pymatsolver.direct import Mumps -from simpeg_drivers.electromagnetics.time_domain import ( - TDEMForwardOptions, - TDEMInversionOptions, -) from simpeg_drivers.electromagnetics.time_domain.driver import ( TDEMForwardDriver, TDEMInversionDriver, ) +from simpeg_drivers.electromagnetics.time_domain.params import ( + TDEMForwardOptions, + TDEMInversionOptions, +) from simpeg_drivers.params import ActiveCellsOptions from simpeg_drivers.utils.testing import check_target, setup_inversion_workspace from simpeg_drivers.utils.utils import get_inversion_output