Use LinearImplicitSystem/LinearSolver with Nedelec_One variable type

[qianhanhua]

Hi, recently I tried to build model to calculate EM scattering problems using Libmesh. I used the Nedelec_One variable type, which resulted say x dofs and y number of nodes. I defined the system as Linear implicit system as I tried to solve the problem directly by assembling the system matrix K and rhs F and solving the matrix equation Kx=F. The results seem to be wrong and as I checked, I found the size of the solution vector to be 3*y, instead of x as I would expect.

Would like to ask gurus here, if it is right to use linear implicit system to solve such problem with Nedelec_One variable type? Thanks a lot!

[pbauman]

LinearImplicitSystem should be just fine. I suggest having a look at vector_fe_ex3 (2D) and vector_fe_ex4 (3D); these solve a simple curl-curl equation, but should be similar enough for you to get going, I think.

[qianhanhua]

LinearImplicitSystem should be just fine. I suggest having a look at vector_fe_ex3 (2D) and vector_fe_ex4 (3D); these solve a simple curl-curl equation, but should be similar enough for you to get going, I think.

Thank you Paul. I actually followed vector_fe_ex1 to construct my code because I thought it was more straightforward. Since it didn't work, I will try to follow vector_fe_ex3 and ex4 as my next step.

[jwpeterson]

It would probably be most helpful if you could post or link to your code, but one comment is I think that NEDELEC_ONE has DOFs at mid-edge nodes but not vertices (?), so you have to use one of the second-order element types for your problem. For example, when I run vector_fe_ex4 I get these numbers:

TET10
n_nodes()=2585
n_elem()=1536
n_dofs()=2156

HEX20
n_nodes()=425
n_elem()=64
n_dofs()=300

HEX27
n_nodes()=729
n_elem()=64
n_dofs()=300

Note that the HEX20 and HEX27 have the same number of DOFs. Is there any chance that you are calling all_second_order() to convert your Mesh, and the number of nodes is not what you were expecting?

[qianhanhua]

Hi, John, thank you for the prompt reply. My code looks similar to vec_fe_ex1, except for variable type that I assigned to NEDELEC_ONE. For mesh I used TET10 which should be fine.

As you said, the dof for NEDELEC_ONE resides on the edges, so I expect to have the same number of unknowns as the number of dofs, instead of number of the nodes. e.g. The final system matrix K that was assembled is having a size say 2156 * 2156 (using the example you mentioned). So I expect the solution vector to have a size of 3 * 2156 for 3D problem, but what I obtained was 3 * 2585 which corresponds to the number of nodes. That lead to my question as first mentioned.

So to me there are three possibilities:

1. To solve for NEDELEC_ONE variable type, I should not use LinearImplicitSystem which assumes nodal unknowns. I should switch to the way vec_fe_ex3 and ex4 were doing.
2. LinearImplicitSystem does not assume nodal unknowns, there are some settings that I didn't set properly so that libmesh didn't set the size of solution vector right.
3. It is right to get the solution size 3*2585, because it was the interpolated result after the system got solved. Error of my code resides in the matrix assembly process (though I thought the assemble function was fine after several round of debugging).
   Could you comment which possibility is more possible?

Below attaches how I defined the equation system and obtained size of the solution vector:

equation_systems.add_system<LinearImplicitSystem> ("Scatter");
equation_systems.get_system("Scatter").add_variable("u", FIRST, NEDELEC_ONE);
equation_systems.get_system("Scatter").attach_assemble_function (assemble_scatter);
equation_systems.init();
equation_systems.print_info();
equation_systems.get_system("Scatter").solve();
std::vector<Number> soln;
equation_systems.build_solution_vector(soln);
out<<"solution size is "<<soln.size();

ps: I created my mesh using gmsh, read directly into libmesh without any all_second_order() conversion.

[jwpeterson]

OK, we need to be a bit careful about the distinction between different kinds of "solution" vectors. The System solution vector,

equation_systems.get_system("Scatter").solution

will be System::n_dofs() in length, i.e. matching the numbers I posted in my previous comment. When you build the "nodal" solution vector, as you are doing in the code snippet above, in 3D that will be of length 3*mesh.n_nodes() as you observed. The distinction between the two is that the former is used for the actual linear algebra (assembly/solve) aspects of the code while the latter is used for visualization (e.g. in Paraview). The latter is formed by evaluating the former at the mesh nodes and using the finite element basis function expansions.

[qianhanhua]

Thank you John! That indeed cleared my doubts about the solution vector size. So is it possible to use the LinearImplicitSystem to deal with Nedelec_One variable types, or must we follow what was presented in vector_fe_ex4?

[jwpeterson]

As @pbauman mentioned, it is fine to use LinearImplicitSystem for your code. vector_fe_ex4 uses "FEMSystem"-style assembly, which is newer than (and in many ways preferable to) the older style LinearImplicitSystem assembly code, but either is capable of assembling the NEDELEC_ONE FEs.

[qianhanhua]

@jwpeterson Thank you for the quick answer. Thank @pbauman too. You guys made my life less miserable. Now I am certain that the error with my existing code should reside in the assembly function. Let me jump back to that and at meanwhile give vector_fe_ex4 a try.

[qianhanhua]

Hi, sorry I am here again. After a full day debugging I am running out of leads, thus I want to ask for some further help here. I was trying to build a FEM model for mie scattering using libmesh. So basically it is about a y-polarized plane wave traveling in vacuum along x axis and scattered by a dielectric sphere with diameter 1/3 of the wavelength. The computation domain was truncated by a PML surrounding box.

The governing equation is Curl(mu_inverse*Curl(E_scatter)) - k_0*k_0*epsilon_r*E_scatter = k_0*k_0*(epsilon_r-1)*E_inc, where epsilon_r is the relative permittivity of the sphere. Bellow attaches the essential sections of the assemble function, other unattached preparation steps are similar if not identical to vector_fe_ex1.

for (const auto & elem : mesh.active_local_element_ptr_range())
{
     dof_map.dof_indices (elem, dof_indices);
     fe->reinit (elem);
     Ke.resize (dof_indices.size(), dof_indices.size());
     Fe.resize (dof_indices.size());

     int domain_id = elem->subdomain_id();			
     if (domain_id == 3)
     {
	 epsilon = calc_PML_tensor(*elem, Lpml);
         mu_inv = epsilon.inverse();
     }
     for (unsigned int qp=0; qp<qrule.n_points(); qp++)
     {	
        /***** Fill in Ke matrix w.r.t. governing equation*****/
        for (std::size_t i=0; i<phi.size(); i++)
            for (std::size_t j=0; j<phi.size(); j++)
	    {
		if (domain_id == 1)
			Ke(i,j) += (curl_phi[i][qp]*curl_phi[j][qp]- k_0*k_0*phi[i][qp]*epsilon_r*phi[j][qp])*JxW[qp];
		else if (domain_id == 2)
			Ke(i,j) += (curl_phi[i][qp]*curl_phi[j][qp]- k_0*k_0*phi[i][qp]*phi[j][qp])*JxW[qp];
		else		  		
             		Ke(i,j) += (curl_phi[i][qp]*mu_inv*curl_phi[j][qp]- k_0*k_0*phi[i][qp]*epsilon*phi[j][qp])*JxW[qp];
	     }
        /***** Fill in right hand side, only elements within the sphere have non-zero RHS *****/
	Number phase(cos(k_0*q_point[qp](0)), -sin(k_0*q_point[qp](0)));

        if(domain_id==1)
	{
	     for (std::size_t i=0; i<phi.size(); i++)
                       Fe(i) += JxW[qp]*k_0*k_0*(epsilon_r-1)*Ay*phase*phi[i][qp];
	} 
     }
     /***** Impose boundary condition, u=0 for all boundary edges *****/
     {
	 ElemSideBuilder side_builder;
        for (auto s : elem->side_index_range())
           if (elem->neighbor_ptr(s) == nullptr)
           {
              const Elem & side = side_builder(*elem, s);
              dof_map.dof_indices(&side, dof_indices_side);

	      for (auto idx_side: index_range(dof_indices_side))
		    for (auto idx_elem: index_range(dof_indices))
		    {
			if (dof_indices[idx_elem]==dof_indices_side[idx_side])
		        {
				const Real penalty = 1.e10;
				Ke(idx_elem, idx_elem) += penalty;
                         }
	            }
            }
       }
       matrix.add_matrix         (Ke, dof_indices);
       system.rhs->add_vector    (Fe, dof_indices);
    }
}

I wonder whether there is any misuse of libmesh library or obvious errors (which my eyes were too blind to see) with this assembly procedure so that when I call
equation_systems.get_system("Scatter").solve(); the program is able to finish running but produces wrong results (When plotted in visit, the pseudo-color plot of magnitude of u does not look like a typical mie scattering pattern). Much appreciate if you can take a look!

[pbauman]

You'll need a good preconditioner. I suggest starting with a direct solver to ensure you've implemented your equations correctly and then you can look into the literature for better preconditioners. With PETSc, you can do -ksp_type preonly -pc_type lu and if you installed a package like superlu with PETSc, then also -pc_factor_mat_solver_type superlu to enable a direct solver.

[qianhanhua]

Thank you Paul! Right now I want to buy a ticket to TX and buy you guys a lunch. Now the results looks more or less correct (at least pattern looks familiar). I will go read about preconditioners.

[qianhanhua]

Hi, Paul and John, after several rounds of refining the FEM simulation, the pattern of the plotted results now looks good and I want to verify the result against analytical solutions. I am thinking about getting the difference between FEM results and the analytical solution at each node, then output the difference results to exodus file format so that I can plot the difference distribution over the entire mesh.

Is it possible for me to deduct the analytical solution from the solution vectors directly? Or is there a better way to implement this?

[pbauman]

You can use the ExactSolution class. See for example https://github.com/libMesh/libmesh/blob/devel/examples/vector_fe/vector_fe_ex4/vector_fe_ex4.C#L141

[qianhanhua]

Hi, Paul, thank you for the reply. I found that ExactSolution class only reports several metrics for errors, e.g. L2-Norm error (Let me know if I am wrong). But I wanted to visualize the distribution of errors. So I came out with the following code:

std::vector<Number> soln, soln1;
equation_systems.build_solution_vector(soln);

/* Function to calculate exact solution at point (x,y,z)*/
ExactValue exact_sol;
VectorValue<Number> Eacc, Delta;
int idx=0;

for(auto node: mesh.node_ptr_range())
{	
        Point pt = *static_cast<Point *>(node);
	Eacc = exact_sol(pt(0), pt(1), pt(2));
        VectorValue<Number> Efem(soln[3*idx], soln[3*idx+1], soln[3*idx+2]);
        Delta = Efem-Eacc;
	soln1.push_back(Delta(0));
	soln1.push_back(Delta(1));
	soln1.push_back(Delta(2));
        idx++;
}
std::vector<std::string> names;
equation_systems.build_variable_names (names, nullptr, nullptr);
ExodusII_IO(mesh).write_nodal_data("delta.e", soln1, names);

This way I am able to output the difference between FEM solution and analytical solution into a ExodusII file and check its distribution. Just wanted to share here for any comments and for use of others who may need it.