London Equation bugfix

CHANGELOG.md

@@ -45,6 +45,8 @@ The format of this changelog is based on
   - Added adaptive time-stepping capability for transient simulations. The new ODE integrators
     rely on the SUNDIALS library and can be specified by setting the
     `config["Solver"]["Transient"]["Type"]` option to `"CVODE"` or `"ARKODE"`.
+  - Fix bug in London equation implementation where a curl-curl term was added to the
+    stiffness operator instead of a mass term.
 
 ## [0.13.0] - 2024-05-20
 

docs/src/config/solver.md

@@ -428,10 +428,14 @@ iterative solver choices, and the default choice depends on the iterative solver
   - `"Default"`
 
 `"DivFreeTol" [1.0e-12]` :  Relative tolerance for divergence-free cleaning used in the
-eigenmode simulation type.
+eigenmode simulation type. Ignored if non-zero
+[`config["Domains"]["Materials"]["LondonDepth"]`](domains.md##domains%5B%22Materials%22%5D%5B%22LondonDepth%22%5D)
+is detected.
 
 `"DivFreeMaxIts" [1000]` :  Maximum number of iterations for divergence-free cleaning use in
-the eigenmode simulation type.
+the eigenmode simulation type. Ignored if non-zero
+[`config["Domains"]["Materials"]["LondonDepth"]`](domains.md##domains%5B%22Materials%22%5D%5B%22LondonDepth%22%5D)
+is detected.
 
 `"EstimatorTol" [1.0e-6]` :  Relative tolerance for flux projection used in the
 error estimate calculation.

palace/drivers/eigensolver.cpp

@@ -154,7 +154,7 @@ EigenSolver::Solve(const std::vector<std::unique_ptr<Mesh>> &mesh) const
   // Construct a divergence-free projector so the eigenvalue solve is performed in the space
   // orthogonal to the zero eigenvalues of the stiffness matrix.
   std::unique_ptr<DivFreeSolver<ComplexVector>> divfree;
-  if (iodata.solver.linear.divfree_max_it > 0)
+  if (iodata.solver.linear.divfree_max_it > 0 && !space_op.GetMaterialOp().HasLondonDepth())
   {
     Mpi::Print(" Configuring divergence-free projection\n");
     constexpr int divfree_verbose = 0;

palace/models/spaceoperator.cpp

@@ -821,7 +821,7 @@ void SpaceOperator::AddStiffnessCoefficients(double coeff, MaterialPropertyCoeff
   // Contribution for London superconductors.
   if (mat_op.HasLondonDepth())
   {
-    df.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetInvLondonDepth(), coeff);
+    f.AddCoefficient(mat_op.GetAttributeToMaterial(), mat_op.GetInvLondonDepth(), coeff);
   }
 }