Help with Defining Bilinear Form using NGSolve

Main Forum Old-Forum - do not post here
1

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:

  1. Does my code for the weak form correctly replicate the variational problem I have? My code works but I wanted to double-check 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

  2. 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

2

Hi,
Latex was buggy in this forum last time, thats why we disabled it. I just checked and it seems to be fixed so I enabled it again.

This are just some remarks, not a complete guidance :wink:

First of all your problem is not symmetric, if you want to solve a nonsymmetric problem you shouldn’t declare the matrix symmetric, you should use either umfpack or pardiso (not sparsecholesky) for any matrix inverse and not the CG as an iterative solver. For example GMRes works in the non symmetric case as well.

Second, the trace operator of the HCurl space only gives you the tangential component (only the continuous part), so from your F only the tangential component will be taken. Further if you specify dirichlet boundaries, these dofs will be eliminated from the system. So specifying everything as dirichlet and then putting in boundary terms on the rhs doesn’t make sense.

What do you mean with your code works? Best is usually setting up some problem with an analytic solution and comparing it. You can as well set up a mesh with only very few elements and compute expected element matrices and vectors by hand and compare it to the ones you get.

The small L2 term in the demo is used to get rid of the kernel of the operator (all gradients). Does your operator have a kernel?

I’m not sure how you derived your equations and what your goal is. Is it even well posed? Your equations are not coupled so you could just solve the second one and use the solution in the first one. The second one is symmetric and could be solved similar to the tutorials. Then from the first you get something like

\int_\Omega \nabla \times B B_t dx = \int_\Omega (jcA + Q_2) B_t dx - \int_{\partial \Omega} F_2 B_t dS

For this to be solveable there have to be some restrictions on Q_2 and F_2…

Hope this helps :slight_smile:

Best
Christopher

3

UPDATE: My apologies, but there was a typo in the LaTeX equations I provided in my first post: in the second equation, the first term should be j (dB) . A_t instead of j (dA) . A_t.

Dear Christopher (and anyone else),

Thank you for your response! It was more than helpful in developing my understanding! I apologize for the delay in my response - I was moving to a different city; but am back with my studies.

If you (or anyone else) has the time, I would sincerely appreciate thoughts on the following questions I had about your guidelines:

i) Am I still allowed to simultaneously solve the two equations using NGSolve? Or would it be best to solve the second one first and plug it into the first as you said? I do not understand why the two equations are not coupled. Let’s say I multiply the first equation by -1 and have the constant c absorb the sign, then the curl term will be positive. Would this make an NGSolve solution easier - and will GMRes still make the most sense?

ii) I found only one online document about GMRes (in the tutorial: 5.5 FETI-DP in NGSolve IV: Inexact FETI-DP). The demo calls DPSpace_Inverse, which takes in sparsecholesky, but I am not sure how to incorporate either umfpack or pardiso. Unless, these two should be tested out without a solver (i.e. a.mat.Inverse(fes.FreeDofs(), inverse=“pardiso”) * f.vec). Regardless, I figured I should replace the CGSolver code in my code example by

GMRESSolver(mat=a.mat, pre=c.mat, precision=1e-6, printrates=False),

but I have two issues with this:
a) Is it an issue that I am working with complex numbers, but GMRESSolver does not accept complex = True, unlike CGSolver? What is a good alternative for nonsymmetric systems that involve complex numbers, otherwise?
b) Do you have any hints on choosing which preconditioner to use (local, direct, multigrid, h1amg, bddc, or a custom one), except by trial and error? I was not able to figure one out for this system of equations on paper and was hoping NGSolve had a workaround. I tried running the solver without a preconditioner, but received an error.

iii) In the LaTeX equations you provided, is the cross product of B or B[sub]t[/sub]? I was just wondering how you derived it, especially what you used to switch the order of the two terms in the product (integration by parts again?).

iv) I did not ask this earlier, but I am wondering whether you had advice on how to properly treat the equations if the coefficients c and d were 3 x 3 constant tensors instead (motivated by material coefficients). Is it trivial? When I try defining them as matrix CoefficientFunctions, I get a dimensional mismatch error in SymbolicBFI. How would I go about tackling this with NGSolve? I tried something similar using the Python package FEniCS, and it seemed as if a possible solution was to iterate over each matrix element.

v) To answer the questions you had asked:

  • By “code works”, I simply meant that it ran and gave me numerical solutions; without any regard to precision or accuracy.
  • I do not know whether this problem is well-posed. My main focus is trying to understand NGSolve better. But, I got these equations by considering a very broad form of Maxwell’s equations from electromagnetism. I am hoping to test NGSolve out with some research problems, but want to understand its capabilities better.
  • As for whether my operator has a kernel, I am not sure because I haven’t theoretically looked at these equations too much. I’m not sure whether there’s a way to numerically check this by setting source terms = 0.

Thank you again for all the care and time you put into your previous reply. They settled several doubts I had for a while!