Greetings! I have been testing and using the PGI formulation in SimPEG as the basis for my master’s thesis and I have some questions I was hoping to get answers to.

The Solver problem:

I was using an earlier version of SimPEG which when used to run the inverse problem gave me the expected results. The output for the inversion code said **Done using same solver and solver_opts as the Simulation3DIntegral problem.** But now, for some reason, the Simulation3DIntegral class has its solver and solver opts set to None, which results in the default Pardiso solver being used for my problem. This in turn changes the outputs to my inversion for the exact same data. How may I go about rectifying this issue?

The problem of misfits

This is related to the previous question in a way. The geophysical data misfits at each iteration seem to be converging well for the gravity data, but plateau after a certain point for the magnetic data. This happens for an inversion I had successfully run before, with the exact same p.f. data, uncertainties, petrophysical data and their statistical distributions. I don’t know why this is happening, but I have the results from each of my tests documented and the difference is clear- the magnetic model is not being recovered. Do you see any reason that might be happening?

The problem of cell-by-cell lithological information prior

I have been trying to incorporate geological information as a prior in the WeightedGaussianMixture class, but so far I have been unsuccessful. Do I have to use a different function to make the model or is there something I am missing?

Solver: You should not have to set a solver for potential fields. Those PF examples in the documentation should show you everything you need with regards to that.

If you have issues balancing your various data misfits, as in the PF examples above, consider using the data misfit reweighting strategy. It is an heuristic strategy to iteratively balance the dataset at each iteration

# iteratively balance the scaling of the data misfits
# initialize values
scaling_init = directives.ScalingMultipleDataMisfits_ByEig(chi0_ratio=[1.0, 100.0])
# iterative reweighting
scale_schedule = directives.JointScalingSchedule(verbose=True)

To include geology through different GMM proportions (gmm.weights_) at various locations (i.e. a geology model), you can pass a numpy array of size (“number of active cells”, “numbers of units”) to gmm.weights_ (like gmm.weights_ = your_array). As they are probabilities of each presence of each rock at each location, each row of your array needs to some to 1, with purely positive values (if you want to set a probability of zero, use something like 1e-128 instead; computations are done in log-space and a zero would trip the computation)

Thank you for the quick reply. I have been using the PF scripts you have linked. I have modelled my scripts on them. The reweighting of data misfits through the chi factor line is also set in my script. It gave results earlier, but for some reason is not working right now. It might be an issue from my side, I will check again. Thank you for the help on the geological information, just to confirm; will the row indices of the array correspond to the spatial indices in the block model?

Sorry to trouble you again, but the problem persists. The misfits are being rescaled at every iteration to give more weight to the magnetic data misfit, but the plateau is not going away, resulting in the model not converging. I’m using the same data and alpha parameters as I used earlier. Following are a few iterations of the previous inversion (before the update);

And the following are some iterations (after the gravity misfit reaches its target) from the inversion I am carrying out now. The data, alpha values and all parameters (bounds, tolerances, etc) are equal in both cases.

The final value of the magnetic data misfit stays around 320000, and the smallness misfit blows up. Do you have any idea as to what might be causing this? Please let me know. Thank you for being so helpful; I appreciate it.