Applying boundary conditions for Maxwell's equations using VectorTools::project_boundary_values_curl_conforming

[Ross Kynch]

Hi,

I'm working on a Maxwell solver, so am using the FE_Nedelec elements. I'm trying to use VectorTools::project_boundary_values_curl_conforming to handle the boundary values, but I'm getting compilation errors when I use the following:

VectorTools::project_boundary_values_curl_conforming(dof_handler, 0, DirichletBoundaryValues<dim>(), 0,

    hanging_node_constraints, boundary_values);

The errors goes away if I remove the final argument, but I'm fairly sure the boundary_value mapping should be included, and the solver (using PETSc) complains about missing diagonal entries if I do. What am I doing wrong here?

How should I be declaring the boundary_values mapping for this call? I'm currently using:

std::map<types::global_dof_index,double> boundary_values;

as one would for the standard finite elements.

Thanks for any help

[Timo Heister]

> VectorTools::project_boundary_values_curl_conforming(dof_handler, 0,
> DirichletBoundaryValues<dim>(), 0,
> hanging_node_constraints, boundary_values);

Take a close look at the documentation and the arguments of that
function. The last argument should be an optional mapping. From the
name of the variable it looks like you are hoping to fill a std::map
with the boundary values, but this is not what the function does: the
boundary conditions will be filled into the ConstraintMatrix instead.


--
Timo Heister
http://www.math.clemson.edu/~heister/

[Ross Kynch]

Thanks for the reply. If I understand correctly, the BCs are put into the ConstraintMatrix, hanging_node_constraints, via project_boundary_values_curl_conforming and so don't need to use MatrixTools::apply_boundary_values as I would for standard elements?

I've adjusted my code so the optional argument is removed, and I don't make a called to apply_boundary_values, but I now get an error from the PETSc solver.. so something is still off. Are there any examples of the use of FE_Nedelec elements available? - I've only found snippets of code from people asking for help, but not code which is shown to be working afterwards.

[Timo Heister]

> Thanks for the reply. If I understand correctly, the BCs are put into the
> ConstraintMatrix, hanging_node_constraints, via
> project_boundary_values_curl_conforming and so don't need to use
> MatrixTools::apply_boundary_values

correct.

> as I would for standard elements?

You can use ConstraintMatrix instead of the std::map approach with all
kinds of elements (and it is the method of choice). See step-6.

> I've adjusted my code so the optional argument is removed, and I don't make
> a called to apply_boundary_values, but I now get an error from the PETSc
> solver.. so something is still off.

Are you using ConstraintMatrix::distribute_local_to_global in your
assembly (again see step-6)? Is your matrix singular?

> Are there any examples of the use of
> FE_Nedelec elements available?

No tutorial step-xy, sadly. @Markus: do you have anything like that
lying around that you might want to share?

[Ross Kynch]

Apologies if I've replied twice.. my first attempt appears to have disappeared.

 

>Are you using ConstraintMatrix::distribute_local_to_global in your 
>assembly (again see step-6)? Is your matrix singular? 

Yes, I use that to assemble the global matrix. I've checked the local matrices and they're only singular with elements of order 0, which I think is to be expected?

 I think I've fixed most of my issues by increasing the quadrature order. The PETSc solver is converging now, although I've not checked the error convergence yet.

>No tutorial step-xy, sadly. @Markus: do you have anything like that 
>lying around that you might want to share? 

That'd be very helpful if there are any. I'll add my example code once I've finished, hopefully it'll be useful to someone else.

Thanks a lot for the help. 

[Markus Bürg]

>> Are there any examples of the use of
>> FE_Nedelec elements available?
>
> No tutorial step-xy, sadly. @Markus: do you have anything like that
> lying around that you might want to share?

Unfortunately, I do not comment my own codes and, thus, have nothing close to being valid for a tutorial right now. However, feel free to assign a tutorial number to it and I will try to converge to a full tutorial over time.

Best,
Markus

[Markus Bürg]

Yes, I use that to assemble the global matrix. I've checked the local matrices and they're only singular with elements of order 0, which I think is to be expected?

No, also Nédélec elements of order zero are perfectly fine to use. They are piecewise constant along edges, but linear into the other tangential directions of the cell.

Best,

Markus

[Ross Kynch]

No, also Nédélec elements of order zero are perfectly fine to use. They are piecewise constant along edges, but linear into the other tangential directions of the cell.

Apologies, I meant the local matrix without any BCs applied. It seems to be behaving as expected now.

One additional query. I'm getting some odd behaviour in the number of CG iterations to converge when running on different numbers of processors using PETSc, although the norm of the error is the same (at least up to 6 or so decimal places - I'd expect some small differences). For example, on 1 processor it is taking ~780 iterations to reach convergence, and on 4 it is taking ~1400 iterations - I've not changed the solver control at all.

Is this normal, or does it point to an error in my code?

Thanks again.

[Timo Heister]

> One additional query. I'm getting some odd behaviour in the number of CG
> iterations to converge when running on different numbers of processors using

> Is this normal, or does it point to an error in my code?

With different processor counts the ordering of the unknowns changes.
This means slight variations are normal.

> For example, on 1 processor it is taking ~780 iterations to reach convergence, and on 4 it is taking ~1400 iterations - I've not changed the solver control at all.

That difference looks too big (as do the absolute numbers, btw.). What
preconditioner are you using? Something like SOR will depend heavily
on the number of unknowns.

General advice: get your program working with deal.II SparseMatrix and
a direct solver like UMFPACK first. Then worry about parallelization
and linear solvers.

[Ross Kynch]

Using the using the BlockJacobi preconditioner and using a tolerance of 1e-12*norm(rhs), which is why the iteration count is fairly high. I've just tried it with the preconditioner disabled and the iteration counts are closer for different numbers of processors (and also a lot lower), so I'll assume that it's down to the ordering changes.

Shall take your advice on the direct solver and parallelisation - perhaps getting a little overexcited in not having to code everything from scratch as I previously worked with spectral element methods for the Navier-Stokes equations and wrote everything myself in F90.

Thanks again for the help and quick replies!

[Wolfgang Bangerth]

On 03/04/2014 12:20 PM, Ross Kynch wrote:
> Using the using the BlockJacobi preconditioner and using a tolerance of
> 1e-12*norm(rhs), which is why the iteration count is fairly high. I've
> just tried it with the preconditioner disabled and the iteration counts
> are closer for different numbers of processors (and also a lot lower),
> so I'll assume that it's down to the ordering changes.

BlockJacobi means that PETSc is doing an ILU decomposition of the
diagonal block of the part each processor "owns" of the matrix. If you
have only a single MPI process in your program, this computes a global
ILU decomposition of the matrix and the preconditioner will be pretty
good. On the other extreme, if you had as many MPI processes as degrees
of freedom, then each process would own a single row, its diagonal block
would be 1x1 and it will only compute the inverse of this one block --
which means that it will do a Jacobi preconditioner, which is known to
be pretty useless.

So, with BlockJacobi, you have to expect that the preconditioner
degrades as the number of processors becomes larger.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Ross Kynch]

Thanks for the information on the preconditioners - that explains all the behaviour I've observed so far.

I have another issue now - hope it's ok to continue to ask here:

I've been considering non-homogenous Dirichlet BCs for the curl-curl equation and am having some odd behaviour with mesh refinement. It appears that every global mesh refinement I perform results in a doubling of the solution value on the boundary. The only thing I can think of is that there is some edge length scaling being applied in the boundary conditions. Is there an additional function call I should be making to account for this?

Apologies if I'm completely wrong, but I've looked through my code and can't find the culprit which leads me to the thinking it is in the application of the boundary conditions.

[Markus Bürg]

Dear Ross,

no, it is support to work by just using VectorTools::project_boundary_values_curl_conforming. Could you please create a small test case to illustrate what your issue is?

Best,

Markus

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Ross Kynch]

Apologies for the delay - was out riding my bike a lot over the weekend!

I've trimmed it down to a simpler problem and removed all the MPI/PETSc stuff. Originally I was solving the benchmark of an L-shape domain using a solution with fractional order Bessel functions to create a singularity at the corner, which was giving this strange doubling of the boundary value in the solution. The problem seems to still occur using a toy solution: I'm getting issues with the solution at the boundary changing with mesh refinement.

So I've basically got a square domain [-1,1]^2 with

curl(curl(u)) + u = f

w/

 u(0) = cos(pi*x)*sin(pi*y) + C

 u(1) = -sin(pi*x)*cos(pi(y) + C

 

 f(0) = (2*pi^2 + 1)*cos(pi*x)*sin(pi*x) + C

 f(1) = -(2*pi^2 + 1)*sin(pi*x)*cos(pi(y) + C

If C=0, then this is fine and I get the expected solution and convergence, etc. However, if C is set to any other value, then the boundary values seem to grow as the mesh is refined - see the values on the edges of the square from the VTK files, each cycle of refinement makes the value grow.

Maybe my experience with Maxwell's equations is not enough, but I don't think this is an invalid solution, and I wouldn't expect this doubling even if I were using the wrong BCs/RHS - surely the solution should still converge to something?

Anyway, see the attached code, I've left C as 0.1 (change bc_constant in the RightHandSide and ExactSolution classes) for now, although maybe it's more obvious that things are doubling with refinement with C=1

Thanks again for the help.

R

[Markus Bürg]

Dear Ross,

please check out the latest subversion. I believe that it should work now.

Best,

Markus

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

<testBC.cc>

[Wolfgang Bangerth]

> So I've basically got a square domain [-1,1]^2 with
>
> curl(curl(u)) + u = f
>
> w/
> u(0) = cos(pi*x)*sin(pi*y) + C
> u(1) = -sin(pi*x)*cos(pi(y) + C
> f(0) = (2*pi^2 + 1)*cos(pi*x)*sin(pi*x) + C
> f(1) = -(2*pi^2 + 1)*sin(pi*x)*cos(pi(y) + C
>
> If C=0, then this is fine and I get the expected solution and convergence,
> etc. However, if C is set to any other value, then the boundary values seem to
> grow as the mesh is refined - see the values on the edges of the square from
> the VTK files, each cycle of refinement makes the value grow.
>
> Maybe my experience with Maxwell's equations is not enough, but I don't think
> this is an invalid solution, and I wouldn't expect this doubling even if I
> were using the wrong BCs/RHS - surely the solution should still converge to
> something?
>
> Anyway, see the attached code, I've left C as 0.1 (change bc_constant in the
> RightHandSide and ExactSolution classes) for now, although maybe it's more
> obvious that things are doubling with refinement with C=1

Ross -- with Markus's changes, this is the output I now get:

.......................
Cycle 0:
Number of active cells: 16
Number of degrees of freedom: 840
Cycle 1:
Number of active cells: 64
Number of degrees of freedom: 3280
Cycle 2:
Number of active cells: 256
Number of degrees of freedom: 12960
Cycle 3:
Number of active cells: 1024
Number of degrees of freedom: 51520
cycle cells dofs L2 Error
0 16 840 1.32580012e-04
1 64 3280 4.19612400e-06
2 256 12960 1.31545840e-07
3 1024 51520 4.11406879e-09
.........................

Is this correct? Assuming that it is, I've made this into a test for the
regression testsuite. Thanks a lot for the program!

Markus -- thanks a lot for fixing this!

[Ross Kynch]

Apologies for not replying sooner, thought I had! Yes this is all correct and it also worked for an analytical solution with a corner singularity (h/p refinement yielding similar results to my test codes in matlab). Many thanks or fixing the issues so quickly to Markus and I'm glad my efforts are of some use. After some detours in my work, I'm now working on a 3D Maxwell solver for use in inverse problems, so sure I'll be back for more help in the coming months!

Ross

[Che Gl]

Hi Ross,

i'm trying to solve some application with Maxwell equation e.g. curl u^-1 curl u + ku = j or more complex curl u^-1 curl u + ku = j - kappa*(du/dt) and i have looked into your code testBC.cc as posted above.
By assembling the global matrix i'm confusing a bit, as you did:

                    cell_matrix(i,j) += ( fe_views.curl(i,q_point)[0]*fe_views.curl(j,q_point)[0]
                                         + dotprod(value_i,value_j) )*fe_values.JxW(q_point);

But as definition of curl u = (1/det Fk)*Fk* curl u^, inserting this into the integral curl u * curl v * dx --> (1/J) * (Fk curl u^) * (Fk curl v^) * dx^ and in 2D we have (1/J) * (curl u^) * (curl v^) * dx^
Here i'm confusing the way you assembled the global matrix can be right for 2D and 3D due to different factor (with Fk and without Fk) ?
Furthermore the term (1/J) was not considered in the equation, instead  (xJ).

Did you test your code similar to testBC.cc for any application? or could you give a small example as you mentioned in another thread?
I would be really thankful with any explanation and suggestion.

Best,

To

[Ross Kynch]

Hi Che,

That was construction was just done quickly for my 2D test problem. I've recently been considering the complex valued 3D problem and have switched to using extractors. My current matrix assembly is as follows:

if (block_index_i == block_index_j)

  cell_matrix(i,j) += ( (1.0/current_mur)*(  fe_values[E_re].curl(i,q_point)*fe_values[E_re].curl(j,q_point)

                                           + fe_values[E_im].curl(i,q_point)*fe_values[E_im].curl(j,q_point) )

      + current_kappa_re*(  fe_values[E_re].value(i,q_point)*fe_values[E_re].value(j,q_point)

                                           + fe_values[E_im].value(i,q_point)*fe_values[E_im].value(j,q_point) )

     )*fe_values.JxW(q_point);

else

  if (block_index_i == 0)

    cell_matrix(i,j) += -current_kappa_im*(  fe_values[E_re].value(i,q_point)*fe_values[E_re].value(j,q_point)

  + fe_values[E_im].value(i,q_point)*fe_values[E_im].value(j,q_point)

 )*fe_values.JxW(q_point);

  else if (block_index_i == 1)

    cell_matrix(i,j) += current_kappa_im*(  fe_values[E_re].value(i,q_point)*fe_values[E_re].value(j,q_point)

                                          + fe_values[E_im].value(i,q_point)*fe_values[E_im].value(j,q_point)

                                         )*fe_values.JxW(q_point);

Where block_index_i and j correspond to the real or imaginary parts and E_re and E_im are the FEValuesExtractors for the real/imaginary parts of the system. Also note I've split the real/imaginary parts of the (-omega^2*epsilon + i*omega*sigma) parameter which I denote Kappa (it'd be k in your equations).

I'm not quite following your point about the transformation of curl from the global to reference element. I was under the impression that those terms are included in the use of the curl function from FEValuesViews (which I access via the extractor), if I'm wrong then I'd be happy to know!

Hope my code helps you, anyway.

[Wolfgang Bangerth]

> I'm not quite following your point about the transformation of curl from the
> global to reference element. I was under the impression that those terms are
> included in the use of the curl function from FEValuesViews (which I access
> via the extractor), if I'm wrong then I'd be happy to know!

Yes. Whenever you call
fe_values[E].value(i,q)
or
fe_values[E].curl(i,q)
then these are values and curls in *real* space, not reference space. This
does not preclude the possibility of bugs in computing these quantities, but
the intend at least is that you would never have to worry about the
transformations at all.

[Che Gl]

Many thanks for your code Ross and your explanation Wolfgang.
I were not sure about the point of "fe_views.curl(i,q_point)", i thought it just shows the component of curl operation without any transformation and it should look like "fe_values.shape_curl (i, q_point)" similar to " fe_values.shape_grad (i, q_point)" in the example. Anyway now i understand the point with transformation, since all transformations were already done and included in FEValue class.
But hereby i have another question: what is the difference between
+ fe_values.shape_grad (i, q_point) and fe_views.gradient (i, q_point) and furthermore  using something like fe_values[].gradient(i,q_point) with FEValuesExtractors. I have tested the second option for step-3 by using fe_views.gradient (i, q_point) and got the same solution.

+ Furthermore i have another question according to method using to solve curl curl A = j (in 2 D) with A is magnetic vector potential, this equation is actually not the same problem as that one from Ross:  curl([mur]^-1 curl(u)) + u = f, with f is a vector having 2 or 3 components
    - as we are in 2D considering a plane X-Y, there exist only 1 compoment of A and j, mean Az and jz (only ). How can the problem can be treated (as 3D spacedim?), since we have 2D problem but A and j have the third component not in the plane we are working on...
   - What is the best way to compute and visualize e.g. B-Field if Az was already computed, since we have Bx= (dAz/dy) and By= (dAz/dx).

I would be really happy with any suggestion!!!

Best

To

[Wolfgang Bangerth]

> But hereby i have another question: what is the difference between
> + fe_values.shape_grad (i, q_point) and fe_views.gradient (i, q_point) and
> furthermore using something like fe_values[].gradient(i,q_point) with
> FEValuesExtractors. I have tested the second option for step-3 by using
> fe_views.gradient (i, q_point) and got the same solution.

Correct. I think you will understand the difference if you read through the
documentation model "Vector-valued problems".

> + Furthermore i have another question according to method using to solve curl
> curl A = j (in 2 D) with A is magnetic vector potential, this equation is
> actually not the same problem as that one from Ross: curl([mur]^-1 curl(u)) +
> u = f, with f is a vector having 2 or 3 components
> - as we are in 2D considering a plane X-Y, there exist only 1 compoment
> of A and j, mean Az and jz (only ). How can the problem can be treated (as 3D
> spacedim?), since we have 2D problem but A and j have the third component not
> in the plane we are working on...

The 2d and 3d versions of the problem are not the same although we often use
similar notation for them when writing it in mathematical terms. This is the
same as the plane-strain or antiplane-strain formulations of elasticity in 2d.
You will have to treat 2d and 3d problems slightly differently in your case,
as it is not entirely trivial to write them in a dimension-independent way.

> - What is the best way to compute and visualize e.g. B-Field if Az was
> already computed, since we have Bx= (dAz/dy) and By= (dAz/dx).

So B=grad Az. This should be straightforward to visualize using a
DataPostprocessorVector object. You'll have to write a class that is derived
from this class.

[Che Gl]

Hi Wolfgang,

thanks for your advice

Correct. I think you will understand the difference if you read through the
documentation model "Vector-valued problems".

 where could i find that documentation 

>     - What is the best way to compute and visualize e.g. B-Field if Az was
> already computed, since we have Bx= (dAz/dy) and By= (dAz/dx).

So B=grad Az. This should be straightforward to visualize using a
DataPostprocessorVector object. You'll have to write a class that is derived
from this class.

No.   B=curl A with A = (0,0,Az)T,  so B has only 2 components Bx= (dAz/dy) and By= (dAz/dx). What is the best way you would suggest to compute e.g. Bx= (dAz/dy) and By= (dAz/dx).
With FEEvaluationAccess< dim, dofs_per_cell_, n_q_points_, dim, Number >::get_curl?

Best.

To

[Wolfgang Bangerth]

> Correct. I think you will understand the difference if you read through the
> documentation model "Vector-valued problems".
>
> where could i find that documentation

http://www.dealii.org/developer/doxygen/deal.II/group__vector__valued.html

> > - What is the best way to compute and visualize e.g. B-Field if Az was
> > already computed, since we have Bx= (dAz/dy) and By= (dAz/dx).
>
> So B=grad Az. This should be straightforward to visualize using a
> DataPostprocessorVector object. You'll have to write a class that is derived
> from this class.
>
> No. B=curl A with A = (0,0,Az)T, so B has only 2 components Bx= (dAz/dy) and
> By= (dAz/dx). What is the best way you would suggest to compute e.g. Bx=
> (dAz/dy) and By= (dAz/dx).

The same as what I suggest above -- using DataPostprocessor.


> With FEEvaluationAccess
> <http://www.dealii.org/8.1.0/doxygen/deal.II/classFEEvaluationAccess.html><

> dim, dofs_per_cell_, n_q_points_, dim, Number >::get_curl?

No, the FEEvaluationAccess class has a completely different purpose.