Simple simulation of electrical potential / current.

[Drinausaur]

Hello,

I am not a FEM expert, but an electronics guy who reads about FEM online ( including the MFEM guide ).

I am trying to set up a "simple" simulation of some pieces of copper, add voltage sources ( DBC if I understand properly ), and current sources ( NBC ).

I started with the voltage sources starting from example 1, but I can't figure out why it does not work. ( I have been looking at that code for days now ).
It seems that I am able to set the voltage to some regions, but not all, and the regions that are assigned take different values.

Before that code, I load the mesh, store voltage values in m_dbcv and the region they correspond to in m_dbcs

    new( &m_fec ) mfem::H1_FECollection( 1, m_mesh.Dimension() );

    new( &m_fespace ) mfem::FiniteElementSpace( &m_mesh, &m_fec );
    Report( "Number of unknowns: " + std::to_string( m_fespace.GetTrueVSize() ) );

    // Set conductivity
    mfem::ConstantCoefficient conductivity( 1 / 1.68e-8 ); // 1 / ( ohm.m )
    mfem::ConstantCoefficient one(1.0);
    mfem::ConstantCoefficient zero(0);

    // Set the Dirichlet values in the solution vector
    mfem::GridFunction u(& m_fespace );
    u = 0.0;

    Report( "Number of element attributes: " + std::to_string(m_mesh.attributes.Size()) );

    // Prepare the source term in the right-hand-side
    mfem::LinearForm b(&m_fespace);
    b.AddDomainIntegrator(new mfem::DomainLFIntegrator( one ));

    // Define the Dirichlet Boundary conditions
    mfem::Array<int> ess_tdof_list;

    for (int i=0; i< m_dbcs.Size(); i++)
    {
        Report( "Applying Dirichlet on region " + std::to_string( m_dbcs[i] ) + " : " + std::to_string( m_dbcs[i] )  );
        mfem::Array<int> dirichlet_attribute =
        dirichlet_attribute.Append( m_dbcs[i] );
        int dirichlet_value = m_dbcv[i];
        mfem::Array<int> dirichlet_tdof_list;
        m_fespace.GetEssentialTrueDofs( dirichlet_attribute, dirichlet_tdof_list);
        ess_tdof_list.Append( dirichlet_tdof_list );
        mfem::ConstantCoefficient bc_coeff(dirichlet_value);
        u.ProjectBdrCoefficient(bc_coeff, dirichlet_tdof_list);
    }

    b.Assemble();

    // Prepare the bilinear form
    mfem::BilinearForm a(&m_fespace);
    a.AddDomainIntegrator(new mfem::DiffusionIntegrator( conductivity ));
    a.Assemble();

    // Prepare the linear system with enforcement of the essential boundary
    // conditions
    mfem::SparseMatrix A;
    mfem::Vector B, X;
    a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);

    // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner.
    mfem::GSSmoother M(A);
    mfem::PCG(A, M, B, X, 1, 200, 1e-12, 0.0);

    // Recover the solution
    a.RecoverFEMSolution(X, b, u );
    u.Save("sol.gf");

And then I get this: ( I tried to set each element to a different voltage, I did not connected them, yet )

here is my mesh:
pcb_3d.msh.zip

Could you help me understand why it does not work?

[mlstowell]

Hello, @Drinausaur,

Unfortunately, this is a common pitfall when working with attributes. We need to document this better. The array you pass to GetEssentialTrueDofs() is not meant to be the set of attributes of interest. Instead it is a "marker" array. Its length should be mesh.attributes.Max() and it should contain zeros and ones where the ones indicate the attributes of interest by their position in the array.

That being said, there is another issue: GetEssentialTrueDofs() works with "element" attributes as opposed to "boundary" attributes. I'm not sure which you intended. Also, ProjectBdrCoefficient does not expect true DoF indices but rather a "boundary attribute marker array". This is a "marker" array as described above but its length is mesh.bdr_attributes.Max().

I would suggest avoiding GetEssentialTrueDofs() because it probably isn't needed in this context. Start by creating a marker array using either:

mfem::Array<int> dirichlet_marker(m_mesh.attributes.Max());
dirichlet_marker = 0; dirichlet_marker[m_dbcs[i]-1] = 1;
u.ProjectCoefficient(bc_coeff, dirichlet_marker);

or

mfem::Array<int> dirichlet_bdr_marker(m_mesh.bdr_attributes.Max());
dirichlet_bdr_marker = 0; dirichlet_bdr_marker[m_dbcs[i]-1] = 1;
u.ProjectBdrCoefficient(bc_coeff, dirichlet_bdr_marker);

With one of these inside your loop you should set the proper values in u.

Best wishes,
Mark

[Drinausaur]

Thank you for your answer,

I made some changes to my code with your help, and it certainly eliminated some mistakes, but the simulation is not yet as expected. I suppose I messed up with the equations.

    new( &m_fec ) mfem::H1_FECollection( 1, m_mesh.Dimension() );

    new( &m_fespace ) mfem::FiniteElementSpace( &m_mesh, &m_fec );
    Report( "Number of unknowns: " + std::to_string( m_fespace.GetTrueVSize() ) );

    // Set conductivity
    mfem::ConstantCoefficient conductivity( 1 / 1.68e-8 ); // 1 / ( ohm.m )
    mfem::ConstantCoefficient one(1.0);
    mfem::ConstantCoefficient zero(0);
    mfem::GridFunction phi(& m_fespace );
    
    Report( "Running solver... " );
    // Initialize the electric potential with its boundary conditions
    phi = 0.0;
    

    // Prepare the source term in the right-hand-side
    mfem::LinearForm b(&m_fespace);
    b.AddDomainIntegrator(new mfem::DomainLFIntegrator( one ));

    mfem::Array<int> dbc_bdr_attr(m_mesh.bdr_attributes.Max());
    for ( int i=0; i < m_dbcs.Size(); i++ )
    {
        Report( "Trying to set voltage on region " + std::to_string( m_dbcs[i] ) + " : " + std::to_string( m_dbcv[i] ) );
        mfem::ConstantCoefficient voltage(m_dbcv[i]);
        mfem::Array<int> dirichlet_bdr_marker(m_mesh.bdr_attributes.Max());
        dirichlet_bdr_marker = 0; dirichlet_bdr_marker[m_dbcs[i]-1] = 1;
        phi.ProjectBdrCoefficient(voltage, dirichlet_bdr_marker);
    }

    b.Assemble();

    // Prepare the bilinear form
    mfem::BilinearForm a(&m_fespace);
    a.AddDomainIntegrator(new mfem::DiffusionIntegrator( conductivity ));
    a.Assemble();

    mfem::Array<int> ess_tdof_list(0);
    mfem::Array<int> ess_bdr(m_mesh.bdr_attributes.Max());
    ess_bdr = 0;

    // Prepare the linear system with enforcement of the essential boundary
    // conditions
    mfem::OperatorPtr A;
    mfem::Vector B, X;
    a.FormLinearSystem(ess_tdof_list, phi, b, A, X, B);

    // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
    mfem::GSSmoother M((mfem::SparseMatrix&)(*A));
    mfem::PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
      
    // Recover the solution
    a.RecoverFEMSolution(X, b, phi );
    phi.Save("sol.gf");

Just for some context, I am working on an open source PCB editor software ( KiCad ) and we'd like to add some FEM capabilities.

For now, the thing we'd like to run is a simple voltage / current simulation, where we can set voltage sources ( so dirichlet BC ) and current sources ( Neumann BC I suppose ) and then let the resistivity parameter determine the gradient electric potential everywhere.

So in that example here, what I'm trying to achieve is setting some "pillars" to a potential ( so they should be uniformly colored ), and the current that results should create a potential gradient along connections.

We have made some successful test with another solver, but we'd like to test mfem.

And this is the equation we'd like to solve: ( this is from my notes when working with the other solver )
At the moment, we focus on DC conduction as this is the certainly the simplest FEM simulation we can run.

[Drinausaur]

Could you suggest a working preconditioner/solver for that specific problem that would work with the serial version?

Going with the parallel version in the future is not excluded, but at the moment, I'd need something that works on basic cases.
Moreover, our FEM interface might need some rework to support mfem more easily, might as well do that soon.

The idea is to fix things on our side before adding things that could break it even further such as going parallel ^^'

[jandrej]

Direct solvers are probably your best bet and are usually very robust until you reach a certain system size. MFEM is integrated with UMFPACK and MUMPS. Both should be easy to install with a package manager and then you can link them after enabling the option in the MFEM build configuration.

[Drinausaur]

sorry, when I wrote solver, I meant those that are integrated in mfem.
The idea is to use mfem in order to limit the number of dependencies. ( therefore using mfem directly )

My question was, after reformulation:

What should I replace those lines with in order to hope for convergence? ( or is there another thing I should replace ? )

    // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
    mfem::GSSmoother M((mfem::SparseMatrix&)(*A));
    mfem::PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);

[jandrej]

Oh I see. You don't want to introduce any TPL. If that is the case your best choice is what you have there. PCG with Gauss-Seidel, which is unfortunately weak for larger Laplace like systems. Increase the number of iterations (here you have 200) to something very large like 5000 and try to figure out a good initial guess. E.g. when you expect your solution to have a minimum potential of 3.0V then you can set your initial guess in the GridFunction to that.

[v-dobrev]

If PCG+Gauss-Seidel starts to get slow as you refine the mesh, you can then try the geometric multigrid solver (assuming you do the mesh refinements in MFEM and not before that in your meshing tool). See examples/ex26.cpp for an example of using the geometric multigrid for a diffusion problem.

[mlstowell]

I ran your mesh in a modified version of ex1p and the results look fine.

While it is an unusual mesh with rather extreme elements I don't think this is the problem. I haven't looked closely at the code yet but perhaps we've missed something there.

[Drinausaur]

I have been doing some tests for the past few days, and the simulation is working.
However, I have some current that leaks out of the simulation.
In the picture below, I set a circle to 1V and another one to 5V. There should be a current between the two ( and therefore potential changes ), but not elsewhere.

I tried setting a Neumann BC to 0, but it is not what I want to enforce. I don't want the net current exiting the simulation to be 0, I want the current exiting the simulation to be 0 everywhere. What is the fem / mfem function I should use for that ?

[jandrej]

If you have a "grounded outer layer", your boundary condition on the surface is Dirichlet = 0V. In 3D, depending on what your model should represent, then this means to also set the surface facing us and the one in the opposite to be zero. That would mean you can only see the potential if you slice into the visualized field.

[mlstowell]

To have zero current flow at a boundary the condition you want is $n\cdot\sigma\nabla\phi=0$. This is a homogeneous Neumann boundary condition and it happens to be the default boundary condition when using H1 finite elements with this operator so you don't need to do anything. I think @jandrej is on the right track, it looks like your outermost surfaces are being set to 0V using a Dirichlet Boundary condition. Remove these boundaries from your Dirichlet BC and it will probably work as expected.

[Drinausaur]

Oh, this was (again) a convergence problem. The code did converge, but the tolerance was not tight enough and lead to this strange result --'. Sorry