Pardiso on Windows: difficulty to deal with the null space

[thast]

In the current version SimPEG.EM.Static.DC, cell-centered formulation, the mix boundary condition does not take care of the null space for Neumann. Solving it on Linux and Mac does not create any issue but the solution obtained on various Windows machine are off by an order of magnitude. Using Joe version of pyMKL and downgrading to MKL 11.3.1 did not affect this behaviour.

link to the script: https://www.dropbox.com/s/e1jjzm08tiwaab1/DCR_Poisson_CC.py?dl=0

image: we see a relative errors of more than 100% (1e0) up to 1000% (1e1) when run on a window machine (relative error is about 1e-1 o mac on Linux, image will be upload later)

More details: On my own machine, the result is wrong first time. On @micmitch machine, he got the same answer as Mac and Linux when run first time but get the wrong answer if the script is rerun...

@lheagy ; @rowanc1 ; @sgkang

[micmitch]

Updated code so that the null space is explicitly dealt with by zeroing the 1st row and A and setting A[0,0] =1. Perhaps that will help ease the discrepancies...

[rowanc1]

Nice.
On Wed, Sep 14, 2016 at 7:23 PM Michael Mitchell notifications@github.com
wrote:

Updated code so that the null space is explicitly dealt with by zeroing
the 1st row and A and setting A[0,0] =1. Perhaps that will help ease the
discrepancies...

—
You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub
#429 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AA3vYWqEZOONWLLQ5l2gLEk5byeJpA4Fks5qqKwLgaJpZM4J8QAj
.

[thast]

Yes @micmitch , that is what I did in my version and it works. My point was more to emphase the difference of result between Mac-Linux and Windows. (and also emphase that if something works on a system, it might not on another one if the null space is not handle)

Also, I did not push this change because it would deactivate all @sgkang 's work on mix boundary conditions.

Last point, I just checked it and the problem is the same on the Master branch! So Windows users would experience this problem using the DC CC formulation from the master branch! My opinion would be to fix the null space at least on Master (A[0,0] += 1) while we figure out a smarter way that does not cancel @sgkang works. What do you think @rowanc1 ?

[micmitch]

Not sure why it cancels @sgkang work with BC. While I was working on it I also put in an option to allow the user to choose whether they wanted Neumann or Dirchlet BC to be used for the CC problem. Testing with analytics indicates that it is working well. This work is all on DC/dev.

[lheagy]

Just had a look at Adam's paper on this. They set A[0,0] += 1/V[0,0]. This might scale the matrix better - may be worth a quick check. (http://library.seg.org/doi/pdf/10.1190/1.2402499?ct=a7309ccf7453e003c35de383a22e8591dc46f73cb9275fe1c34c463c22fa50a5e0437017b4445a5dbb332aae5fe482b08ce36bf8a4512572d48b57ae4b47743f)

[thast]

@lheagy : true. A[0,0] has no reason to be close to 1 and +1 can be a huge perturbation. In finite volume +=1/V[0,0] make more sense indeed

@micmitch: Nice, I will take a look at it then. my concern with mixed boundary conditions was that the matrix perturbation would not make sense if we use a pure Derichlet BC (and I do not know about it when we use a mixed BC...)

[micmitch]

For Problem_CC functionality has been added so that the user can specify Dirchlet or Neumann BC. When using Neumann we take care of the null space as suggested by Lindsey (reference Adam's article above).

i.e. Zero first row of A and set A[0,0] = 1/V[0,0]

From analysis that Thibaut and I have done, we are wondering if we should set Dirchlet instead of Neumann as the default BC for the CC problem. Dirchlet generally produce results which are closer to the analytic solution. Dirchlet even appears to behave nicer when we reduce the number of padding cells substantially.

[sgkang]

In Adam's formulation matrix A do not include volume term, whereas SimPEG's solution does includes volume. So when assigning a value to A[0, 0], we do not need 1/Vol[0].
Mixed boundary condition, which treats this issue is implemented. See #649 #650