Support for 3D PDE

[BBBBBbruce]

Hi,

I am trying to use Zombie to solve a 3D Poisson equation; I wonder if you have any 3D examples? it seems I need to change all the files in demo/ since it is a hard-coded 2D example.

Thank you

[rohan-sawhney]

We don't currently have 3D demos in this repository, but it's high on the priority and the plan is to release some in the coming weeks/months (can't commit to an exact timeframe unfortunately).

However, the library itself fully supports 3D equations, and the structure of the code for 2D vs 3D problems remains pretty much unchanged.

In C++, you can reuse most of the boiler plate code to run the solver in the model_problem.h and demo.cpp files - you'll basically have to change the template DIM parameter from 2 to 3.

The python bindings also support both 2D and 3D problems, and again the demo.py can be pretty much be reused. Since python does not support templates, function calls come with explicit _2d and _3d tags. You can view the full API as follows:

>>> import zombie
>>> help(zombie)

[rohan-sawhney]

I've templated the code in demo.cpp on dimension (and value type). For scalar valued PDEs in 3D, you should be able to call runSolver<float, 3>.

[BBBBBbruce]

Hi Rohan,

thank you, I have read the source code and the main problem for me to transfer to 3D is to define the boundary condition. Currently in 2D the BC is parameterised as a 2D image, which for 3D would be a bit difficult. A similar way is to simply stored as a tensor. But using a tensor is not quite efficient in Eigen. Since zombie.PDE.source/.dirichlet/.robin, etc are functions. I directly changed the lambda functions in model_problem.setupPDE(), where I given the analytical boundary definition. For example,

    pde.dirichlet = [this, bMin, maxLength](const Vector3& x, bool _) -> float {
        float value = 0;
        if (x[0]<0.1 || x[0]> 0.9) {
            value = 1.0f;
        }else{
            value = 0.0f;
        }
        return value;
    };

But this seems not quite a natural way to use this in C++ because I need to recompile everytime I change my input. But using function input would be much natural in python. So, that is what I am currently working on. I could share the progress after I implement this if you think that is useful.

[rohan-sawhney]

Yep, the C++ interface is designed so that users can programmatically define the PDE via lambdas, but the downside is that code needs to be recompiled. Future versions of the library will allow the user to define the PDE in a shader, so that will help bypass recompilation.

I'd be interested in hearing how you plan to programmatically define the PDE in Python. The python bindings currently have callbacks to define PDE functions, e.g., constant or (2D/3D) dense grid based boundary conditions. See here. Is there a more flexible way you have in mind to programmatically define the PDE in Python?

Also note that the library provides a DenseGrid class to query arbitrary data in 2D or 3D, which could be useful for defining PDE callbacks in either dimension. Naturally, such a grid can become memory intensive in 3D, especially if it has a high resolution.

[BBBBBbruce]

The shader does sounds like a good plan, but I am not familiar with coding, so I will leave that to you.

I did noticed in the demo.py you have DenseGrid class, but I think the lambda function the cpp demo is a better inferface.

I found pybind11 can be used to pass function input to c++; a simple example is:

# c++
#include <pybind11/pybind11.h>

namespace py = pybind11;

py::object example(int x, py::function func) {
    py::object result = func(x);
    //bool b = result.cast<bool>();
    return result;
}

PYBIND11_MODULE(my_module, m) {
    m.def("example", &example);
}

c++ -O3 -Wall -shared -std=c++11 -fPIC $(python3 -m pybind11 --includes) my_module.cpp -o my_module$(python3-config --extension-suffix)

# py
import my_module

def my_func(x):
    return x *x-2

result = my_module.example(5, my_func)
print(result)

I only tested this when we build the demo in python because my other codes are also in python.

Additionally, from what I checked, this can also be use for cpp demo, but more tricky. The logic is to load the python code as string and calling the python terminal to execute it. Since the python demo is sufficient for my case, I didnt test this.

[rohan-sawhney]

Yep, it's possible to pass along C++ functions to python and vice versa---the dense grid PDE callbacks are a version of the former in the zombie python codebase using nanobind.

As you mention, it's also possible to pass python functions to C++ using nanobind, however, this does come with a performance overhead. Specifically, the C++ zombie codebase is multithreaded but when calling a python function from C++, nanobind (and pybind) invoke the Global Interpreter Lock despite each thread not needing any locking, which I suspect can cause a slowdown.

[BBBBBbruce]

you are right. performance is indeed what I need to keep for comparisons. I will kept the cpp code for now.

BTW, I notice a problem. the queries.insideDomain(pt) would classify the node on the positive boundary as insideDomain which leads to segment fault.

This is a common corner case for the odd-even rule for raycasting. but I dont know why that would causing segment fault. I tried a unit square from (0,0) to (1,1). any node of (1,x) or (x,1), where 0<x<=1 would lead to segment fault. The code I traced to is:

In zombie/deps/fcpw/include/fcpw/core/primitive.h

  bool findClosestPoint(BoundingSphere<DIM>& s, Interaction<DIM>& i, bool recordNormal=false) const {
      int nodesVisited = 0;
      return this->findClosestPointFromNode(s, i, 0, this->pIndex, nodesVisited, recordNormal);
  }

I suspect the this->findClosestPointFromNode is not initialised, because I cannot step further in debugging.

The reason I would having the boundary is because I am comparing to the FEM solution and therefore the gridpoints are points from traignulated mesh. so, the boundary points are included. My current solution to this is to remove the boundary element first. which adds quite a lot overheads tho.

Do you know how can I handle this error properly?

[rohan-sawhney]

Thanks, I'm surprised calling insideDomain is resulting in a call to findClosestPoint. If you look at the implementation of insideDomain here, there is no call to a closest point function, instead we are performing raycasting.

The solvers (WoS, WoSt, BVC, etc.) never actually use the insideDomain function. Are you getting a segfault from simply calling this function?

For sample points on the boundary, you can explicitly set the type parameter here to either SampleType::OnAbsorbingBoundary or SampleType::OnReflectingBoundary. That should be sufficient as an input to the solvers.

[BBBBBbruce]

Hi, I get the dirichlet and source worked out. the results looks amazing. However, I am having some trouble setting up Neumann. The observation is that the integral of neumann is zero. and during my debug, I found the problem is that my walk is always find a way to escape, except those directly walked on to the dirichlet.

I also found that the walk can escaped, where this is returned:

       if (!state.onReflectingBoundary && queries.outsideBoundingDomain(state.currentPt)) {
            if (walkSettings.printLogs) {
                std::cout << "Walk escaped domain!" << std::endl;
            }

            return WalkCompletionCode::EscapedDomain;
        }

In my pure dirichlet experiment, this can also happen, but result still correct. when I add Neumann, all non-zero totalReflectingBoundaryContribution walk would resulted in escaped domain.

I wonder should I expect walk escape at all? do you have any insight on this problem maybe?

[BBBBBbruce]

hmmm, I have recreate the problem with a simple example.

I use the exact 2D demo. and change geometry.obj to a unit square.obj

the unit square.obj is:

o boundary
v 1.0 1.0
v 0.0 1.0 
v 0.0 0.0 
v 1.0 0.0 
l 1 2 
l 2 3
l 3 4
l 4 1

and toggle

"runSingleThreaded": true,
"printLogs": true

then the terminal is printing many
Walk escaped domain!

[rohan-sawhney]

Hi, if the domain is closed/watertight, and you want to evaluate the solution in the interior of the domain, then walks should not escape.

I think the issue you are facing might have to do with how you are specifying the input geometry in your obj file. Specifically, in 3D, the input domain should be represented as a triangle mesh, but in 2D, the input domain should be a polyline mesh representing the boundary of the domain (and not the domain itself). See for instance the l flag in the geometry.obj used for the demo.

[BBBBBbruce]

okay, thanks. I managed to get all the 2D cases(using walk on star and enable Dirichlet, source and pure neumann) working where I turned off the flip-orientation flag. why this might causing issue? for 2D case it just swap the order of edge indices. Is it flip the normal?

However, for 3D I still couldnt get it to work. For 3D I inputed the water-tight surface mesh stored in the typical .obj file created from blender and it only has vertex coordinates and faces. I tried to use trimesh to fixed normals, turning flip-orientation on and off. but cannot worked it out.

I attached the results, where left is solved with fem and right is wost, I evaluated at all the points from the mesh for comparison. additionally, I turn on solve on both side, so I can evaluate the points on the boundary. (the boundary is defined by the .obj and internal nodes are tetrahedralised by gmsh and load them instead of uniformly sampling)
(BTW, I checked the interior values, not mapped with the fem solution as well, the color is just indication where it is globally normalized.)

[rohan-sawhney]

Can you share code snippets for loading the geometry and setting up the solver? That'll help me more easily figure out the issue with your code setup. Thanks.

[BBBBBbruce]

Hi thanks.

for PDE setup I am using analytical functions:, so it is pure Neumann, and the boundary condition is setted along z axis, where

Dirichlet: z < 0.01
Neumann 0: 0.01 <= z <= 0.99
Neumann 1: z > 0.99

my mesh is normalised in range (0,0,0) - (1,1,1)

    pde.source = [this](const Vector3& x) -> float {
        return 0.0f;
    };

    pde.dirichlet = [this](const Vector3& x, bool) -> float {
        if (x[2]<0.01){
            return (std::sin(2*3.1416*x[0])+1);
        }else{
            return 0.0f;
        }
    };

    //std::cout<<"Debug: ModelProblem::setupPDE: setting up boundary conditions: others"<<std::endl;
    pde.robin = [this, bMin, maxLength](const Vector3& x, const Vector3& n, bool _) -> float {
        if(x[2]>0.99){
            return 1.0f;
        }
        else{
            return 0.0f;
        }
    };
    pde.robinCoeff = [this](const Vector3& x, const Vector3& n, bool _) -> float {
        return 0.0f;
            
    };
    pde.hasReflectingBoundaryConditions = [this, bMin, maxLength](const Vector3& x) -> bool {
        if(x[2]>=0.01){
            return true;
        }
        else{
            return false;
        }
    };

the loading is default.
zombie::loadBoundaryMesh<3>(filename, vertices, segments);
with options:

        "normalizeDomain":false,
        "flipOrientation": false,
        "solveDoubleSided": true

the mesh this: made by meshio(change txt to obj).
pillerNew.txt

the nodes I queried:
wos_eval.csv

using load gridpoint instead of create it uniformly:

void loadGridPoints(std::string gridPath, std::vector<Vector3>& gridPoints)
  
{

    std::ifstream in(gridPath);
    if (!in.is_open()) {
        throw std::runtime_error("Failed to open file: " + gridPath);
    }
    std::string line;

    while (std::getline(in, line)) {
        std::stringstream ss(line);
        std::string field;
        float x, y, z;


        std::getline(ss, field, ','); x = std::stof(field);
        std::getline(ss, field, ','); y = std::stof(field);
        std::getline(ss, field, ','); z = std::stof(field);

        gridPoints.emplace_back(Vector3(x, y, z));
    }

}

[rohan-sawhney]

Thanks, I'll take a closer look at this data and respond soon.