PCG solver did NOT converge in 50000 iterations in a coupled microstrip simulation

[CosimoMV]

I tried to simulate a coupled microstrip line, but no matter what I do, I always get this error: PCG solver did NOT converge.

I tried to simulate the copper interface with PEC and a metal layer with finite conductivity, I varied the parameters "Lc", "DivFreeMaxIts", "EstimatorMaxIt" and "MaxIts" , I tried the solvers "SuperLU", "STRUMPACK" and "MUMPS" with both version 0.12 and 0.13, but I always get this error. I think I am doing something trivial wrong (the system I am simulating is simple), but I do not understand what the problem is.

This is an example of my .json files simulation.json

This is my mesh: tutto.zip

[sebastiangrimberg]

Hi @CosimoMV,

Looking at your mesh, it is likely that your problem could be solved by using a nicer initial mesh with more isotropic elements and smoother size grading between elements. This can be handled in Gmsh by applying size fields, see for example this tutorial. This should help the iterative solvers in Palace to converge.

[CosimoMV]

Hi @sebastiangrimberg , thanks for your help. I tried to recreate the mesh as you said, but I was only getting absurdly large files. I therefore reconstructed the geometry of the coupled microstrip directly with gmesh and meshed it with these parameters:

Mesh.MeshSizeFromPoints=1;
Mesh.MeshSizeFromCurvature=1;
Mesh.MeshSizeExtendFromBoundary=2;
Mesh.ElementOrder=4;

In this way I was able to create files that were quite small and allowed Palace to converge. The problem is that the solution, apart from being unrealistic, is highly dependent on the algorithm chosen to mesh.
In the graph:

are shown the |S21| (dB) vs frequency (GHz) parameters obtained with various meshing algorithms: the first curve name indicates the algorithm for surfaces, the second the one used for volumes, and the third is the name of Palace solver. The S-parameters values, apart from one case, are very different and in any case unrealistic. In all cases, the Palace solver converged.
What do you think is the error? Please tell me if I should close this thread and open another one. By the way, the Palace version I used was v0.12.0-390-g66d4a16.

[sebastiangrimberg]

Hi @CosimoMV, thanks for your help with the mesh generation.

I think I can predict what is going on with your simulation, though I am not 100% certain and cannot test my theory without your configuration file and mesh. There are two issues, I suspect, both related to the use of wave ports in your model:

1. Boundary conditions for the wave port mode solve: When you use a wave port, a boundary mode eigenvalue solve will take place for each wave port. This solve only supports PEC and PMC boundary conditions, and by default all PEC boundaries are inherited from the PEC boundaries of the full 3D model. If you specify a surface impedance or surface conductivity boundary condition for your metal, this won't be automatically interpreted as a PEC for the 2D mode solve. Instead, you need to list these attributes explicitly as "WavePortPEC". See also this PR which proposes at least automatically adding the surface conductivity boundary conditions as PEC by default. Adding this to your configuration file you should should see the correct boundary modes being identified.
2. Mode degeneracy: For the wave ports you have added to your coupled microstrip model, there are two modes, even and odd. The mode which the wave port selects is just given by the "Mode" index value, which selects the desired mode based on the eigenmodes sorted by propagation constant. It could be for some frequencies in a sweep that you will get the even or odd, which will of course make your S21 or S11 plots look strange. One way to enforce this is to split your single wave port into two and enforce even vs. odd symmetry with a PEC or PMC boundary condition on the edge between the ports.

I hope this explanation is helpful. Using wave ports can be complicated and if you see room for improvement in the implementation we would always welcome a PR. It is also worth mentioning that you can visualize the computed wave port modes for each wave port in the ParaView output for the boundary mesh, which should allow you to validate and debug your model and results.

[CosimoMV]

Hi @sebastiangrimberg, first of all thank you very much for your reply and for your help (also in the other thread), it was really useful.

Regarding the first point, I correctly used “WavePortPEC” to set the boundaries for the wave port. If I may say, the PR you mentioned would be a good addition to Palace, because from the point of view of newbie (like me), it isn't immediate to understand that you need to set explicitly PEC boundaries for the waveport, if you doesn't want the default PMC boundary.

Regarding the second point: my concerning was that Palace, with the same configuration, but with mesh from different algorithms, converged to very different solutions that in any case were unrealistic.
So I re-meshed all my geometry with a size field and didn't use gmsh Refinemesh (for some reason with this option the mesh was smaller, but was really slower for Palace to converge). In this way I obtained a quite large mesh file (28.2 MB), but with this mesh Palace converged in a reasonable time and in a few steps to “good” s-parameters from 1GHz to 29Ghz.

After 30 Ghz I obtained very big positive values for abs(S11), exactly like in #267 . So, after your reply, I investigated the problem a little further and got these results:

When abs(S11) became positive the propagating eigen-mode at the boundary was probably an evanescent mode, because the wave vector had a big imaginary component:

@ 28GHz Port 1, mode 1: kₙ = 8.737e+02-6.571e+00i m⁻¹
@ 29GHz Port 1, mode 1: kₙ = 9.049e+02-6.810e+00i m⁻¹
@ 30GHZ Port 1, mode 1: kₙ = 2.069e-01-9.061e+02i m⁻¹ <---
@ 31GHz Port 1, mode 1: kₙ = 2.162e-01-8.903e+02i m⁻¹
@ 32GHz Port 1, mode 1: kₙ = 2.256e-01-8.738e+02i m⁻¹

I recorded the field values at two mirror points in relation to the plane of central symmetry, halfway between the bottom of the miscrostrip and the ground plane:

        --------       --------     Z = 0
        
           p1             p2
           
---------------------------------------------- ground        Z = -h   

f (GHz)	Re{E_z[1]} (V/m)	Re{E_z[2]} (V/m)
2.800000000e+01	+5.374040815e+04	-5.370912572e+04
2.900000000e+01	+5.384068796e+04	-5.382424342e+04
3.000000000e+01	+4.357688567e+03	+4.675744413e+03
3.100000000e+01	+4.652318182e+03	+4.913960087e+03
3.200000000e+01	+4.504950000e+03	+4.745971071e+03

The values for Ez at 28GHz and 29GHz for p1 and p2 have the same modulus and opposite sign, as is right for an odd-mode.
For 30 GHz, 31GHz and 32GHz, Ez has the same sign (so, not an odd-mode) and different modulus (even considering numerical errors, only the first digit is the same), so it would appear that they are not even-modes, and this would confirm my hypothesis that they are evanescent modes.

Then for 30 - 50 GHz I set mode 2 instead of mode 1 and I obtained the “right“ S-parameters – abs(S11) for frequencies higher than 30Ghz continuously connected to the part 1 -30Ghz.
In this case (mode = 2), the situation was a mirror image with inverted frequency range of when mode = 1:

@ 28 Ghz Port 1, mode 2: kₙ = 1.889e-01-9.353e+02i m⁻¹
@ 29 GHz Port 1, mode 2: kₙ = 1.979e-01-9.210e+02i m⁻¹
@ 30 GHz Port 1, mode 2: kₙ = 9.362e+02-7.049e+00i m⁻¹ <-------
@ 31 GHz Port 1, mode 2: kₙ = 9.675e+02-7.288e+00i m⁻¹
@ 32 GHz Port 1, mode 2: kₙ = 9.989e+02-7.527e+00i m⁻¹

I recorded the field values in the usual two mirror points exactly as before:

f (GHz)	Re{E_z[1]} (V/m)	Re{E_z[2]} (V/m)
2.800000000e+01	-5.378851852e+03	-5.045164501e+03
2.900000000e+01	+2.938638665e+03	+3.254652357e+03
3.000000000e+01	+5.383748914e+04	-5.385184002e+04
3.100000000e+01	+5.373062827e+04	-5.379738050e+04
3.200000000e+01	+5.357836749e+04	-5.368964052e+04

Finally, regarding this:

One way to enforce this is to split your single wave port into two and enforce even vs. odd symmetry with a PEC or PMC boundary condition on the edge between the ports.

I had also thought of something similar, but I did not know whether it was possible to superimpose the waveport with a PEC strip on the central symmetry plane. So, if I understand correctly, the solution would be to create two ports per side with appropriate boundary conditions instead of one, and then obtain a s4p file at the end?
Would you mind explaining the procedure in more detail?

[sebastiangrimberg]

Thank you very much for the detailed testing and reporting.

First, I did some thinking about evanescent modes and think I have come up with a solution there. I have added a commit to the PR #264 which improves the wave port boundary mode to exclude evanescent modes correctly. Testing this on the CPW example (as in #267), shows that it fixes the issue we observed above 30 GHz with the S-parameters. You should be able to build with this PR and I am hoping this will resolve many of the issues you identified.

Regarding your last question about resolving even/odd mode ambiguity, yes I think you have the right understanding. This is what I had done in the coupled-CPW example here: https://awslabs.github.io/palace/dev/examples/cpw/. Instead of using one port per side, I split it and use two ports per side and get a 4 dimensional scattering matrix. I should note that I am not 100% certain if this 4-port approach is better or even is a fully correct model for the coupled problem but I think it might be worth trying out.

[CosimoMV]

I did the test with the branch of PR #264 and it worked perfectly even with the coupled microstrip! There is no longer any need to change the mode index at 30GHZ. Tank you very much @sebastiangrimberg !

I thought a bit about how to make sure to only found the odd (or even) modes: I think the best way to do is to exploit symmetry, which would also reduce the size of the problem by a factor of 2 (just for the propagating eigen-mode calculation).
Referring to the picture below (I know, it's really ugly), lets say we want to found a propagating odd-mode at a given frequency. We could split the domain in half respect to the plane of symmetry and calculate the mode only in zone 1 with the appropriate boundary conditions (if we want an odd-mode we must impose a PEC boundary on the symmetry line). Then we could just use the symmetry relations between E and H fields components to find their values also in zone 2. In this way we would be sure to have imposed a predetermined symmetry and would have split the size of the problem in two, although compared to calculating the propagation in the 3D model I do not think it is a major improvement.

By the way, do you have any tips for me regarding mesh creation? To use Palace in the real cases that really interest me, I am forced to start from .step files that come out of a CAD. Do you have any tips on how to do the mesh creation, which seems to me to be the most critical part for the success of the simulation.

[sebastiangrimberg]

That's great about the PR! I will get this reviewed and merged ASAP. And correct about using symmetry for even/odd mode simulation, I agree with you there.

For meshing, I understand this can be a complicated part of the workflow. If you are just importing a STEP file, I think you can still get quite a long way in Gmsh using size fields. When you import a STEP, you should be able to still access the various geometric entities (like the trace domain) and add a size field which controls the mesh size based on distance from a domain or boundary. This is is what I do in the CPW example here and it works quite well. You can read more about that in the Gmsh docs, here: https://gmsh.info/doc/texinfo/gmsh.html#Specifying-mesh-element-sizes. There's also a nice example: https://gmsh.info/doc/texinfo/gmsh.html#t10.

Alternatively, I am not sure what you are using to create your 3D CAD (STEP), but there are many unified CAD + meshing softwares which might simplify things. For example using the OpenCASCADE geometry kernel in Gmsh along with meshing functionality can give you access to a lot of programatic ways to control mesh sizing, avoiding the file import/export of a STEP.

[CosimoMV]

Alternatively, I am not sure what you are using to create your 3D CAD (STEP), but there are many unified CAD + meshing softwares which might simplify things. For example using the OpenCASCADE geometry kernel in Gmsh along with meshing functionality can give you access to a lot of programatic ways to control mesh sizing, avoiding the file import/export of a STEP.

My current goal is to make simulations post-layout, so I have to use the output of our electrical-CAD. The only outputs our CAD produces that are useful here are step or odb++ files. For now I was just opening the step files with freecad to do the finishing touches and importing the final files into gmsh. I also tried meshing with freecad, but the meshing options it provides were too limited.

I was already using OpenCASCADE kernel in gmsh because it is needed for boolean operations. I will do as you say and try to impose a mesh size with fields.

Thank you very much @sebastiangrimberg for your help and advice! I would say that this thread can be closed without any problems. If anything, I'll open a new one if I have further problems, although I sincerely hope not!