Hello,
I am relatively new to NGSolve and have been exploring its functionality using demos, etc. To apply all these methods via a problem of my own, I am attempting to numerically solve for the following variational problem (image attached because I could not figure how to enter LaTeX in here):
Above, I am solving for A and B, with test functions A_t and B_t respectively. j is the imaginary unit. c and d are constants. For i = 1 or 2, Q_i are constant vectors. F_i are constant vectors but acting on the boundary.
Although I made the coupled PDEs above by drawing from the terms in Maxwell’s equations, the setup is unphysical. This is a toy model that I want to play with using different parameters and NGSolve features (to study numerical analysis).
However, because I am new to NGSolve, I am not sure whether my implementation below is correct for the two coupled PDEs. I simply added the two equations above and defined the bilinear form. Essentially, this forum post was made to ask for help in validating the way I implemented this toy model (sanity checks, obvious errors, etc):
a = BilinearForm(fes, symmetric=True, eliminate_internal=True)
a += SymbolicBFI(1j*InnerProduct(c*trial_A, test_B)+ 1j*(InnerProduct(d*trial_A, test_A)))
a += SymbolicBFI(InnerProduct(trial_B, curl(test_B)) + InnerProduct(trial_A, curl(test_A)))
f = LinearForm(fes)
f += SymbolicLFI(InnerProduct(Q_2,test_B) + InnerProduct(Q_1, test_A ))
f += SymbolicLFI(InnerProduct(F_2, test_B.Trace()), BND)
f += SymbolicLFI(InnerProduct(F_1, test_A.Trace()), BND)
I have attached a minimally working example below as well. As you can see, I solve the problem over a mixed H(curl) function space, over a unit cube. The code runs without errors.
Given all this, my questions are:

Does my code for the weak form correctly replicate the variational problem I have? My code works but I wanted to doublecheck I am not abusing any NGSolve functionality, or choosing inefficient/incorrect methods. In particular, I am worried about:
a) summing the two PDEs and defining them together (I did not alter signs in the forms above, etc)
b) using Trace() and BND for F_i (is this the way I define that integral?)
c) using “symmetric=True, eliminate_internal=True” in BilinearForm() (should I add a small L2 term like in some demos? when I did, there was no big difference in solutions)
d) any other subtlety that I missed 
Any advice on my choice of solvers/preconditioners for a problem that looks like this? I could not find a similar problem online, so was not sure which was most valid.
I use the “local” preconditioner for CGSolver because others (i.e. “direct”,“multigrid”,“bddc”): give numerical values that are zero; or large error when I try to solve problems for A and B that I already know analytically; or take too long and run out of memory.
Am I doing anything wrong by sticking to “local”? Or is some other preconditioner/solver more suitable for the problem? How do I go about finding the best solution? Any resource you would recommend, or any obvious trick?
Thank you very much for your time and patience!
https://ngsolve.org/media/kunena/attachments/1276/forumexample.py
Attachment: forumexample.py