Möbius Domain Wall Fermion Solver
The Möbius domain wall object, its variants, and solving with it.
Qlua provides bindings to the MDWF Level III library. Some of what the library implements may not be exposed; this can be extended if needed. The existing bindings are described below.
Creating a MDWF object
The main entity provided is the dwf object. It is implemented only for 4-d lattice and Nc=3. If these conditions are not met, the constructor reports an error and Qlua exits with an error code. There are several variants of the Domain Wall fermions available. The most generic case allows the user to define values of \(b_5\) and \(c_5\). The generic dwf object is created from the following parameters
- The gauge field. As usual, a table of four Color Matrices.
- Boundary conditions. A table of four numbers defining factors picked up by the fermion field when crossing between \(L_i - 1\) and \(0\). As all Lua tables, it is 1-based.
- \(L_s\), the size of the fifth dimension.
- \(M_5\), the 5d mass value.
- \(m_f\), the domain wall mass value.
- \(b_5\), values of b parameter, a table of \(L_s\) real numbers.
- \(c_5\), values of c parameter, a table of \(L_s\) real numbers.
Several special cases of the Domain Wall are provided as a convenience. There are Shamir, Borichi, Chiu, and Möbius variants available. See MDWF documentation for their definitions.
Once the dwf object is created, U could be reused.
One can use multiple dwf objects at the same time. Keep in mind, that each dwf object takes a considerable amount of memory, though.
Here we create a Shamir operator with Dirichlet boundary conditions in the T-direction:
DW_dirichlet = qcd.MDWF.Shamir(U,{1, 1, 1, 0}, 16, 1.8, 0.040, 1.0);
Another example: the standard antiperiodic in time, periodic in space operator:
DW = qcd.MDWF.Shamir(U,{1, 1, 1, -1}, 16, 1.8, 0.040, 1.0);
Applying the Möbius Operator
With the dwf object constructed, one can apply both the Dirac operator and its conjugate.
phi = DW:D(psi); -- the Dirac Operator
eta = DW:Dx(psi); -- the conjugated operator
The Möbius operator is five-dimensional, so its argument is not a single field but a table of Ls DiracFermion fields, one per slice of the fifth dimension; the result is a table of the same shape. A plain DiracFermion or a DiracPropagator is an error.
local psi = {};
for i = 1, Ls do psi[i] = S:gaussian_DiracFermion(); end
local phi = DW:D(psi); -- a table of Ls fields, as is the result
Solving the Dirac Equation
To compute the inverse of the Dirac operator, one has to construct a solver first. Two kinds of solvers are provided.
The standard CG solver is created as follows:
std_solver = DW:solver(1e-22, 10000);
Here the first argument is the required precision. The solver uses the red-black preconditioned normal equation, and stops when \(|Mx-b|^2/|b|^2 < \epsilon\). The second argument is the maximum number of iterations CG could use. If the convergence is not reached within that limit, the solver returns anyway, with the best iterate it has.
A solver returns three values: the solution, the residual it actually reached, and the number of iterations it took. Checking the second of those is the only way to tell a converged solve from one that ran out of iterations, and nothing else will tell you:
local x, eps, iters = std_solver(b);
if eps > 1e-22 then
printf("no convergence: %e after %d iterations\n", eps, iters);
end
The second solver uses mixed single/double approach of the BMW group. You create the mixed solver by saying
mx_solver = DW:mixed_solver(1e-10, 1000, 1e-22, 10000);
Here the first argument is precision to use in the single precision solver. The second argument sets the maximal number of iterations in each single precision inverter. Next is the desired precision of the solver. Argument four sets the total number of iterations the solver may use before bailing out. The mixed solver still produces a double precision solution, but it may run faster than the standard solver. In particular, on the Blue Gene one can get about 30% speedup. On other platforms, the mixed solver could be at a disadvantage. Compare the two before embarking on a production calculation.
When applied to a fermion or a propagator, either solver will produce the solution to the Dirac equation. Here we solve the equation both ways and compute the true error:
std_inv = std_solver(psi);
mx_inv = mx_solver(psi);
b2 = psi:norm2();
std2 = (DW:D(std_inv) - psi):norm2();
mx2 = (DW:D(mx_inv) - psi):norm2();
printf("standard solver error: %.2e\n", std2/b2);
printf("mixed solver error : %.2e\n", mx2/b2);
The :norm2() method exists for every lattice field and computes the square of the norm across the active lattice sites.
Note that with \(c_5 \ne 0\), the four-dimensional source must be converted to five dimensions and multiplied by the operator \(-D_-\), defined in Eq. (2.18) of <http://arxiv.org/abs/1206.5214>. The solution can then be converted back to four dimensions. This is needed to ensure \(\gamma_5\)-hermiticity of the four-dimensional solution.
Memory management
The Dirac operators take considerable amount of memory. Though Qlua does proper memory management and reclaims this memory when the operator object is garbage collected, one may want to clean up earlier. The operator object has the close() method that frees all memory used by the operator object and marks the object against future use.
DW:close();
After closing the operator object, one must not call any of its methods or solvers constructed from it. If one does, Qlua performs the usual catastrophic error shutdown.
|
