CURRENT STATUS: What is wrong with this code?

[jackhamel16]

This issue references the Development branch.

The code currently solves the full LLG Equation using a Jacobian Free Newton Krylov (JFNK) solver with GMRES. I recommend reading a short report I wrote explaining what the code does including some results before reading the rest as it probably won't make sense. Note that some of the notation is a little sloppy, but it is mostly sound. You may see mention of them, but Neumann Iterations do not work due to their dissipative nature so only focus on GMRES looking forward. Additionally, the small timescale issue was fixed by using the proper gyromagnetic ratio. The issue now lies with the speed GMRES solves for the next iteration of V (V is the difference between M between JFNK iterations). It takes a veryyyyyy long time. I have done tests with one particle and a purely DC incident field and it takes hours to run for a ~10 microsecond simulation. If you push the timestep size too small you will start to get wrong results as well. I believe the dt minimum is around 10^-13 seconds. Not great.

So what is wrong? Well, I think the issue lies in the jacobian matrix-vector product approximation being used in conjunction with having fields dependent on M (e.g. the self-interaction term). I have checked analytically if the jacobian matvec approximation matches the explicit jacobian dotted with a vector and it does so this leads me to believe there is an issue buried in the code. I tried removing all fields dependent on M or V and the solver worked very fast and was accurate. All this meant was getting rid of the self-interaction, which I don't think should change the final result because in the LLG equation you'll notice that any field aligned with M has zero contribution due to the cross products. I could be wrong about this when we are dealing with the fields due to V. My hypothesis is that a bug resides in the code handling these operations.

[jackhamel16]

Lastly, I am available to try and help whoever is working on this project with questions they may have! Please reach out to me via email: jackhaml16@gmail.com

[jackhamel16]

This will give a brief history of recent tests and developments with this code starting from January 2019.

1) Set incident field to be a time-varying pulse (I think I used a gaussian) and solved the full LLG Equation for one particle using both Euler and PC Integrators. Results were compared to python numerical integration and the results matched well (I don't recall exact errors and don't want to dig for them right now).

2) Changed how M(Now) is solved for in issue 18

3) Next step was to solve the LLG equation for two particles simultaneously. Using the Newton-Raphson method, this failed (blew up still). I went and double checked input parameters first. I realized the magnetization magnitude was actually too low. At this point we took a step back and evaluated the physics, comparing our code with Vitaliy Lomakin's work. This is when the JFNK + Neumann iterations or GMRES was introduced. See issue 19

5) JFNK solver was implemented using Nuemann iterations (not GMRES yet). JFNK solver was tested against an analytical solution. I think this passed, but in the development branch, the JFNK_test unit test is turned off so I think it should be tested again. This shouldnt be hard to do, if I remember correctly.

6) Blow up issues still resided under the right conditions. This was due to using Neumann Iterations. These were replaced with GMRES. GMRES was tested independent of the code and is correct. It also takes Ax as a product but not explicitly the matrix and vector.

7) At this point, my notes are bad/lost on what was done during the month of May 2019. From memory, I discovered that the gyromagnetic ratio was previously wrong, fixing this helped run at larger timesteps because the actual value is smaller than what was used before thus decreasing the M rate of change. I explicitly formed the jacobian of the full LLG and used it in the JFNK solver. The results of this test should be compared to the implicit jacobian again. Lastly, I checked the jacobian approximation to make sure it reduced to the actual jacobian. I am pretty sure it does, but I kinda want to look into this again because I think this gave me insight into what I think the bug in the code is related to.

WHAT TO DO: So where to go from here? First, get the JFNK solver unit test working again. this will be a good starting point and should refresh my memory on how it works. Second, compare the results of the explicit jacobian to implicit jacobian. Last, check the jacobian approximation analytically that it matches the actual jacobian.

[jackhamel16]

Also backup all of your research files, you daredevil

[jackhamel16]

The JFNK unit test would require actual rework of the code or some tricky/advanced unit testing skills that I do not have so I actually abandoned it. I should truncate the llg equation and compare the results to the analytical solution of it to see if the solver works.

The explicit and implicit jacobian methods match with at most ~1e-5 error. I checked on pen and paper if the two methods match and they don't. I originally thought they did but I can clearly see an error term proportional to the fields due to the v vector (difference between magnetization between iterations of JFNK solver). Additionally, there is no speed up gained by using one method over the other so this begs the question: why use the implicit jacobian matvec...