Nodal Poisson with EB geometry

[Alex-iono]

Hi,
I am an amrex newbie so this might be a stupid question. However, I am trying to solve a nodal Poisson problem on a slightly fancier geometry than in the example code.
I want to simulate two charged electrodes instead of a point charge. I managed to generate the geometry using the EB::SphereIF and assigned sigma(high in spheres, low outside), potential(non zero in, zero outside) and charge(same as pot) based on it.
When I try to solve this I get unphysical results. Am I missing something? Maybe someone has a similar snippet flying around...
Cheers,
Alex

[Alex-iono]

poisson.cpp

Sourcefile for context

[WeiqunZhang]

assigned sigma(high in spheres, low outside)

If you use EB, MLEBNodeFDLaplacian assumes the solid object is a perfect conductor with sigma being infinity.

When I try to solve this I get unphysical results.

What result are you getting? What result are you expecting?

[Alex-iono]

Hi,
Good to know the thing with implicit sigma infinity.
This is the final potential solution that I am getting:

For me it looks like the solver basically did nothing. It is the initial guess for the nodal potentials with a slight blurring due to the fact that I display it cell centers.

Yesterday I also played around with the setOversetMask() option for my NodalLaplacian, trying to set the sphere potentials as known DoFs but I either did it wrong or it doesnt help.

Regards,
Alex

[WeiqunZhang]

The source code you provided is insufficient to determine what's going on. For example,

Real reltol = m_interface_ptr.safe_get<Real>("solver.rel_tol");
...
mlmg.solve({&m_phi}, {&m_q}, reltol, 0.0);

Maybe reltol is 0 or even negative? We could not know. So I can only guess.

Maybe you can set verbose=2 and show us the output on the screen?

I would like to help you, but I would appreciate if you can provide a small reproducer. I would not have the time if the code is say 10,000 lines. So if you have a complete reproducer that is small, I will try to look into it.

[Alex-iono]

Hello again,

this is the verbose output:

_
AMReX (25.10) initialized
Finished generating geometry
MLMG: # of AMR levels: 1
# of MG levels on the coarsest AMR level: 1
MLMG: Initial rhs = 1
MLMG: Initial residual (resid0) = 4.778666667e+17
MLCGSolver_BiCGStab: Initial error (error0) = 1
MLCGSolver_BiCGStab: Final: Iteration 7 rel. err. 6.41275517e-05
MLMG: Iteration 1 Fine resid/resid0 = 2.838797316e-05
MLCGSolver_BiCGStab: Initial error (error0) = 6.41275517e-05
MLCGSolver_BiCGStab: Final: Iteration 201 rel. err. 0.007027425829
MLCGSolver_BiCGStab:: failed to converge!!
MLMG: Bottom solve failed.
MLCGSolver_CG: Initial error (error0) : 1.356566611e+13
MLCGSolver_cg: Final Iteration 6 rel. err. 1.879794754e-05
MLMG: Iteration 2 Fine resid/resid0 = 5.336356641e-10
MLMG: Final Iter. 2 resid, resid/resid0 = 255006696, 5.336356641e-10
MLMG: Timers: Solve = 3.936324827 Iter = 3.900357895 Bottom = 3.884867115
AMReX (25.10) finalized
_

Seems to have convergence issues.... This poisson.cpp is basically the entire project source code at the moment.
I will create a repository with the entire project code in a second, so you can rebuild it.

Regards,
Alex

[Alex-iono]

Thats the minimal reproducer:

https://github.com/Alex-iono/amrex_testsim.git

[WeiqunZhang]

diff --git a/src/poisson.cpp b/src/poisson.cpp
index 8aa9e69..c3c0cc4 100644
--- a/src/poisson.cpp
+++ b/src/poisson.cpp
@@ -81,24 +81,7 @@ void AmrexPoisson::generate(void){
 
         // set initial simga distribution
         GpuArray<Real,3> dx = m_geometry.CellSizeArray();
-        for(MFIter mfi(m_sigma); mfi.isValid(); ++mfi){
-
-            const Box& v_bx = mfi.validbox();
-            const Array4<Real>& mf_array_sigma = m_sigma.array(mfi);
-            // assign cell centered sigmas
-            ParallelFor(v_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k){
-                
-                const RealArray center_pos({(i+0.5) * dx[0], (j+0.5) * dx[1], (k+0.5) * dx[2]});
-                // assign cell centered sigma
-                if (sphere1(center_pos)>0){          // in sphere1
-                    mf_array_sigma(i,j,k) = 1e10;
-                } else if (sphere2(center_pos)>0) {  // in sphere2
-                    mf_array_sigma(i,j,k) = 1e10;
-                } else {                                // outside spheres eps0
-                    mf_array_sigma(i,j,k) = 8.85e-12;
-                }
-            });
-        }
+       m_sigma.setVal(8.85e-12);
         // set initial q and phi distributions
         for(MFIter mfi(m_q); mfi.isValid(); ++mfi){
 
@@ -154,7 +137,6 @@ void AmrexPoisson::solve(){
                         {AMREX_D_DECL(LinOpBCType::Dirichlet,
                                     LinOpBCType::Dirichlet,
                                     LinOpBCType::Dirichlet)});
-    linop.setOversetMask(0, m_mask);
     //for (int ilev = 0; ilev <= max_level; ++ilev) {
     linop.setSigma(0, m_sigma);
     //}

This appears to fix your issue. The way you are doing sets sigma to a huge value in cut cells. Note that you don't need to worry about the value of sigma in covered cells at all. It's effectively infinity. You also don't need to set overset mask. That's for a different purpose.

Other comments. It might be more efficient for big runs if you add more multigrid levels.

         EB2::SphereIF sphere2(domain_x/12, {domain_x/2,domain_y/2,domain_z-domain_z/4}, false);
         auto two_spheres = EB2::makeUnion(sphere1, sphere2);
         auto gshop = EB2::makeShop(two_spheres);
-        EB2::Build(gshop, m_geometry, 0, 0);
+        EB2::Build(gshop, m_geometry, 0, 30);
         std::cout << " Finished generating geometry " << std::endl;
 
         // Factory

and

     bool agglomeration = false;
     bool consolidation = false;
     bool semicoarsening = false;
-    int max_coarsening_level = 0;
+    int max_coarsening_level = 30;
     int max_semicoarsening_level = 0;
     int verbose = 2;
     Real reltol = 1e-6; //m_interface_ptr.safe_get<Real>("solver.rel_tol");

[WeiqunZhang]

Looks like this