Added Px and Py components

[jackieyao0114]

This PR should be reviewed and merged after PR#14.

This PR adds Px and Py components by including the corresponding terms in the Landau free energy and the gradient energy. The derivations are shown in the following document. The final terms added in the code are highlighted in rentangle frames in the document. PxPy_PR16.pdf

This PR followed the following steps as listed in Issue #10, repeated here:

1. Initialize a three component MultiFab for P
2. Modify the RHS of Poisson solve to read div(P) - rho instead of 'dP/dz'
3. Modify the RHS of TDGL equation to account for all three components
4. Update time integrator loop to advance Px, Py, and Pz
5. Pass new material parameters for Landau free energy term and gradient energy terms
6. Add diagnostics for all three components.

We can use the following paper as a reference: https://aip.scitation.org/doi/pdf/10.1063/1.1377855

However, it should be emphasized that

1. The diagnostic multifab Plt does not include Px and Py yet, as we want to do a direct multifab comparison to the previous Pz-only simulations
2. The DoubleDPDz and DPDz functions still remain the previous 1st order one-sided derivative for P_BC_flag_hi == 1. As I have noticed well-above machine precision discrepancy with the 2nd order one-sided derivative compared to the current development branch.
3. In ComputeEfromPhi function, the calculation of Ex, Ey and Ez is still hard coded. We should update them with the new kernels in DerivativeAlgorithm.H.

[ajnonaka]

What do you mean by comment 3? ComputePoissonRHS doesn't need to compute E. If you meant ComputeEfromPhi, note that you can't re-use the DPDz routines since they are specially-coded for the P boundary and interface conditions.

[ajnonaka]

Is your intention for this new feature to work for AMREX_SPACEDIM=2? If so, there need to be some ifdefs added and careful use of AMREX_SPACEDIM as opposed to hard-coded "3".

[jackieyao0114]

What do you mean by comment 3? ComputePoissonRHS doesn't need to compute E. If you meant ComputeEfromPhi, note that you can't re-use the DPDz routines since they are specially-coded for the P boundary and interface conditions.

Nice catch. I meant ComputeEfromPhi.

And nope, we are not using any DPDz for E calculation. There are dphidz in the kernel, which we should use for calculating E.

[ajnonaka]

OK; note that "phi" and "E" have different physical boundary conditions, so you can't use the exact same kernel near the boundary.

[jackieyao0114]

OK; note that "phi" and "E" have different physical boundary conditions, so you can't use the exact same kernel near the boundary.

I'd be curious to learn more.

[ajnonaka]

Now that I think about it more, DphiDz is probably correct. DphiDZ IS Ez...

[ajnonaka]

Yeah, we are specifying the phi boundary condition, and since phi and E are related, E should be picking up the correct boundary condition.

[jackieyao0114]

Yeah, we are specifying the phi boundary condition, and since phi and E are related, E should be picking up the correct boundary condition.

Agreed!

[prkkumar]

I tested this PR with GL_RHS_x(i,j,k) = 0.0; GL_RHS_y(i,j,k) = 0.0; in Source/TotalEnergyDensity.cpp, which turns off the evolution of Px and Py terms. I am able to get good agreement with current development branch.

Simulation of a MFIM stack with the following input (input.txt) was performed and the results at time steps 0, 500, and 1000 were compared with the results obtained using development using ./DiffSameDomainRefined3d.gnu.ex infile1=plt00001000 reffile=./PxPy/plt00001000 (this is for timestep = 1000).

At timestep = 0

Level  L2 norm of Error in Each Component
-----------------------------------------------
 Level: 0   
P : 0 
Phi : 0 
PoissonRHS : 0 
Ex : 0 
Ey : 0 
Ez : 0 
holes : 0 
electrons : 0 
charge : 0 

At timestep = 500

Level  L2 norm of Error in Each Component
-----------------------------------------------
 Level: 0   
P : 3.809748327e-42 
Phi : 1.829106044e-41 
PoissonRHS : 6.417883173e-25 
Ex : 8.520991998e-24 
Ey : 0 
Ez : 3.57696661e-24 
holes : 0 
electrons : 0 
charge : 0 

At timestep = 1000

Level  L2 norm of Error in Each Component
-----------------------------------------------
 Level: 0   
P : 2.643108031e-42 
Phi : 1.593229779e-41 
PoissonRHS : 4.985739282e-25 
Ex : 6.255809832e-24 
Ey : 0 
Ez : 3.172988422e-24 
holes : 0 
electrons : 0 
charge : 0 

These differences are very small and I think this verifies the implementation.

Note that this test was done for prob_type = 1, to make sure that initial states are identical. For prob_type = 2, initial states are obtained using a random number generator. We could use a fixed seed in either case in this case, but, I believe this is sufficient.

[jackieyao0114]

This is nice! Thank you for testing it. We can merge it once Andy approves the PR. And we can fix the minor issues (listed at the end of the description ) with future PRs.