NVIDIA / modulus-sym

Framework providing pythonic APIs, algorithms and utilities to be used with Modulus core to physics inform model training as well as higher level abstraction for domain experts
https://developer.nvidia.com/modulus
Apache License 2.0
137 stars 55 forks source link

🐛[BUG]: total surface area can be wrong with stl tessellation samples #92

Open Karl-JT opened 7 months ago

Karl-JT commented 7 months ago

Version

1.3.0

On which installation method(s) does this occur?

Docker

Describe the issue

When using an stl file (sphere for example) instead of primitive geometry, the total area of all sampled points are different from the actual total area. This leads to errors on computing drag or lift when summing up the pressure on the surface. The reason is in tessellation geometry sampling. The area of samples are calculated based on each triangle, but some triangle will have zero samples.

Minimum reproducible example

@modulus.sym.main(config_path="conf", config_name="config")
def run(cfg: ModulusConfig) -> None:
    point_path = to_absolute_path("./stl_files")
    sphere_geo = Tessellation.from_stl(
        point_path + "/sphere.stl", airtight=True
    )

    ns = NavierStokes(nu=1.0, dim=3, time=False)

    # determine inputs outputs of the network
    input_keys = [Key("x"), Key("y"), Key("z")]
    output_keys = [Key("u"), Key("v"), Key("w"), Key("p")]

    flow_net = FullyConnectedArch(
        input_keys=input_keys,
        output_keys=output_keys,
        adaptive_activations=cfg.custom.adaptive_activations,
    )

    nodes = (
        ns.make_nodes()
        + [flow_net.make_node(name="flow_network")]
    )

    domain = Domain()

    interior = PointwiseInteriorConstraint(
        nodes=nodes,
        geometry=sphere_geo,
        outvar={"continuity": 0, "momentum_x": 0, "momentum_y": 0, "momentum_z": 0},
        batch_size=1000,
        lambda_weighting={
            "continuity": Symbol("sdf"),
            "momentum_x": Symbol("sdf"),
            "momentum_y": Symbol("sdf"),
            "momentum_z": Symbol("sdf"),
        },
    )
    domain.add_constraint(interior, "interior")

    force_p = PointwiseMonitor(
        sphere_geo.sample_boundary(100),
        output_names=["p"],
        metrics={
            "total_area": lambda var: torch.sum(var["area"]),
        },
        nodes=nodes,
    )
    domain.add_monitor(force_p)

    # make solver
    slv = Solver(cfg, domain)

    # start solver
    slv.solve()

if __name__ == "__main__":
    run()

Relevant log output

RuntimeWarning: divide by zero encountered in divide
  np.full(x.shape, triangle_areas[index] / x.shape[0])

The monitored total surface is 9866.225, but the actual surface of the sphere is 31415.93.

Environment details

No response

Other/Misc.

No response

Karl-JT commented 7 months ago

One potential solution could be change the following

invar["area"].append(
    np.full(x.shape, triangle_areas[index] / x.shape[0])
)

to

invar["area"].append(
    np.full(x.shape, np.sum(triangle_areas) / nr_points)
)

or to enforce at least one sample for each triangle.