diff --git a/README.md b/README.md index 7a3de73f..17fe3ade 100644 --- a/README.md +++ b/README.md @@ -67,6 +67,22 @@ Alternatively, if you want to skip pysdf installation, you can run the following docker build -t modulus-sym:deploy -f Dockerfile --target no-pysdf . ``` + +### Windows Container +For building windows just use +``` +docker build -t modulus - < Windows.Dockerfile +``` + +### Additional notes + +New modification suppose to work with pip to (you'll just need to install two additional libralies) +``` +pip install trimesh +pip install mesh_to_sdf +``` +Tesselation will work because instead of using pysdf I used mesh_to_sdf (I have not yet implemented sdf derivatives yet, so some functions will not work on windows (It's only a few functions for visualization differences with real data (as far as I know))) + ## Contributing For guidance on making a contribution to Modulus, see the [contributing guidelines](https://github.com/NVIDIA/modulus-sym/blob/main/CONTRIBUTING.md). diff --git a/Windows.Dockerfile b/Windows.Dockerfile new file mode 100644 index 00000000..7a09cb5c --- /dev/null +++ b/Windows.Dockerfile @@ -0,0 +1,11 @@ +FROM nvcr.io/nvidia/modulus/modulus:23.05 + +RUN rm -rf \ + /usr/lib/x86_64-linux-gnu/libcuda.so* \ + /usr/lib/x86_64-linux-gnu/libnvcuvid.so* \ + /usr/lib/x86_64-linux-gnu/libnvidia-*.so* \ + /usr/lib/firmware \ + /usr/local/cuda/compat/lib + +RUN pip install trimesh \ + mesh_to_sdf \ No newline at end of file diff --git a/modulus/sym/eq/pdes/navier_stokes.py b/modulus/sym/eq/pdes/navier_stokes.py index 2bb1f44e..4547f37c 100644 --- a/modulus/sym/eq/pdes/navier_stokes.py +++ b/modulus/sym/eq/pdes/navier_stokes.py @@ -339,6 +339,91 @@ def __init__(self, T, dim=3, time=True): normal_x * T.diff(x) + normal_y * T.diff(y) + normal_z * T.diff(z) ) +class NavierStokesIncompressible(PDE): + '''Custom implementation of incompressible NS equations''' + + name = "navier stokes incompressible" + + def __init__(self, nu, rho=1, dim=3, time=False): + + '''Parameters + ========== + nu : float, Sympy Symbol/Expr, str + The kinematic viscosity. If `nu` is a str then it is + converted to Sympy Function of form `nu(x,y,z,t)`. + If `nu` is a Sympy Symbol or Expression then this + is substituted into the equation. This allows for + variable viscosity. + rho : float, int + The density of the fluid. Default is 1. + dim : int + Dimension of the Navier Stokes (2 or 3). Default is 3. + time : bool + If time-dependent equations or not. Default is False.''' + + self.dim = dim + self.time = time + + # coordinates + x, y, z = Symbol("x"), Symbol("y"), Symbol("z") + + # time + t = Symbol("t") + + input_variables = {"x": x, "y": y, "z": z, "t": t} + if self.dim == 2: + input_variables.pop("z") + if not self.time: + input_variables.pop("t") + + + # kinematic viscosity + if isinstance(nu, str): + nu = Function(nu)(*input_variables) + elif isinstance(nu, (float, int)): + nu = Number(nu) + + # density + if isinstance(rho, str): + raise Exception("rho must be number") + elif isinstance(rho, (float, int)): + rho = Number(rho) + + + u = Function("u")(*input_variables) + v = Function("v")(*input_variables) + if self.dim == 3: + w = Function("w")(*input_variables) + else: + w = Number(0) + + + + #pressure + p = Function("p")(*input_variables) + + # set equations + self.equations = {} + self.equations["continuity"] = ( + u.diff(x) + v.diff(y) + w.diff(z) + ) + + self.equations["momentum_x"] = ( + u.diff(t) + u * u.diff(x) + v * u.diff(y) + w * u.diff(z) + 1 / rho * p.diff(x) + - nu * (u.diff(x).diff(x) + u.diff(y).diff(y) + u.diff(z).diff(z)) + ) + + self.equations["momentum_y"] = ( + v.diff(t) +u * v.diff(x) + v * v.diff(y) + w * v.diff(z) + 1 / rho * p.diff(y) + - nu * (v.diff(x).diff(x) + v.diff(y).diff(y) + v.diff(z).diff(z)) + ) + self.equations["momentum_z"] = ( + w.diff(t) +u * w.diff(x) + v * w.diff(y) + w * w.diff(z) + 1 / rho * p.diff(z) + - nu * (w.diff(x).diff(x) + w.diff(y).diff(y) + w.diff(z).diff(z)) + ) + + if self.dim == 2: + self.equations.pop("momentum_z") class Curl(PDE): """ @@ -507,3 +592,84 @@ def __init__(self, T="T", D="D", rho=1, vec=["u", "v", "w"]): self.equations[str(T) + "_flux"] += ( Symbol(v) * n * rho * T - rho * D * n * g ) + + +class Euler(PDE): + '''Custom implementation of compressible Euler equations''' + def __init__(self, rho=1, dim=3, ratio_heats=1.4, time=False): + ''' + Parameters + ========== + rho : float, Sympy Symbol/Expr, str + The density of the fluid. If `rho` is a str then it is + converted to Sympy Function of form 'rho(x,y,z,t)'. + If 'rho' is a Sympy Symbol or Expression then this + is substituted into the equation to allow for + compressible Navier Stokes. Default is 1. + dim : int + Dimension of the Navier Stokes (2 or 3). Default is 3. + time : bool + If time-dependent equations or not. Default is False. + ratio_heats : float + Ratio of specific heats. Default is 1.4 (for air). + ''' + self.dim = dim + self.time = time + + # coordinates + x, y, z = Symbol("x"), Symbol("y"), Symbol("z") + + # time + t = Symbol("t") + + # make input variables + input_variables = {"x": x, "y": y, "z": z, "t": t} + if self.dim == 2: + input_variables.pop("z") + if not self.time: + input_variables.pop("t") + + u = Function("u")(*input_variables) + v = Function("v")(*input_variables) + if self.dim == 3: + w = Function("w")(*input_variables) + else: + w = Number(0) + + + # pressure + p = Function("p")(*input_variables) + + # density + if isinstance(rho, str): + rho = Function(rho)(*input_variables) + elif isinstance(rho, (float, int)): + rho = Number(rho) + + + self.equations = {} + + e = p /(rho * (ratio_heats - 1)) + # Energy + E = rho * (e + (u **2 + v **2 + w **2) / 2) + + self.equations["continuity"] = ( + rho.diff(t) + (rho * u).diff(x) + (rho * v).diff(y) + (rho * w).diff(z) + ) + + self.equations["momentum_x"] = ( + (rho * u).diff(t) + (rho * u **2 + p).diff(x) + (rho * u * v).diff(y) + (rho * u * w).diff(z) + + ) + + self.equations["momentum_y"] = ( + (rho * v).diff(t) + (rho * u * v).diff(x) + (rho * v **2 + p).diff(y) + (rho * v * w).diff(z) + ) + + self.equations["momentum_z"] = ( + (rho * w).diff(t) + (rho * u * w).diff(x) + (rho * v * w).diff(y) + (rho * w **2 + p).diff(z) + ) + + self.equations["energy"] = ( + E.diff(t) + (u * (E + p)).diff(x) + (v * (E + p)).diff(y) + (w * (E + p)).diff(z) + ) diff --git a/modulus/sym/geometry/tessellation.py b/modulus/sym/geometry/tessellation.py index ddb14f97..3e9d7e0c 100644 --- a/modulus/sym/geometry/tessellation.py +++ b/modulus/sym/geometry/tessellation.py @@ -20,14 +20,21 @@ import csv from stl import mesh as np_mesh from sympy import Symbol +import trimesh +import threading + + +from mesh_to_sdf import mesh_to_sdf + +#from pysdf import SDF +from copy import copy + + try: import pysdf.sdf as pysdf except: - print( - "Error importing pysdf. Make sure 'libsdf.so' is in LD_LIBRARY_PATH and pysdf is installed" - ) - raise + print("sdf derivatives are not implemented yet for mesh_to_sdf") from .geometry import Geometry from .parameterization import Parameterization, Bounds, Parameter @@ -35,6 +42,65 @@ from modulus.sym.constants import diff_str +def calc_sdf_test(store_triangles, points): + global worked_thread + worked_thread = False + try: + sdf_field, sdf_derivative = libc.signed_distance_field( + store_triangles, points, include_hit_points=True + ) + except Exception as e: + print(f'pysdf doesn\'t work becayse {e}') + + else: + worked_thread = True + + + +def test_pysdf(triangles, invar, airtight): + '''test if PySDF is working and if not uses mesh_to_sdf''' + + if not airtight: + return + + + points = np.stack([invar["x"], invar["y"], invar["z"]], axis=1) + + # normalize triangles and points + store_triangles = np.array(triangles, dtype=np.float64) + minx, maxx, miny, maxy, minz, maxz = _find_mins_maxs(points) + max_dis = max(max((maxx - minx), (maxy - miny)), (maxz - minz)) + minx, maxx, miny, maxy, minz, maxz = _find_mins_maxs(points) + max_dis = max(max((maxx - minx), (maxy - miny)), (maxz - minz)) + store_triangles = np.array(triangles, dtype=np.float64) + store_triangles[:, :, 0] -= minx + store_triangles[:, :, 1] -= miny + store_triangles[:, :, 2] -= minz + store_triangles *= 1 / max_dis + store_triangles = store_triangles.flatten() + + + points[:, 0] -= minx + points[:, 1] -= miny + points[:, 2] -= minz + points *= 1 / max_dis + points = points.astype(np.float64).flatten() + + + + + # Create a separate thread to cath C++ execption + thread = threading.Thread(target=calc_sdf_test, args=(store_triangles, points)) + + # Start the separate thread + thread.start() + + # Wait for the separate thread to complete + thread.join() + + return worked_thread + + class Tessellation(Geometry): """ Constructive Tessellation Module that allows sampling on surface and interior @@ -52,6 +118,15 @@ class Tessellation(Geometry): def __init__(self, mesh, airtight=True, parameterization=Parameterization()): + + #test if can import PySdf + try: + import pysdf.sdf as pysdf + except: + self.pysdf_avalible = False + else: + self.pysdf_avalible = True + # make curves def _sample(mesh): def sample( @@ -111,12 +186,17 @@ def sample( invar["normal_z"] = np.concatenate(invar["normal_z"], axis=0) invar["area"] = np.concatenate(invar["area"], axis=0) + #test pysdf avalibilyty and if not use + self.pysdf_avalible = test_pysdf(mesh.vectors, invar, airtight) + # sample from the param ranges params = parameterization.sample(nr_points, quasirandom=quasirandom) return invar, params return sample + + curves = [Curve(_sample(mesh), dims=3, parameterization=parameterization)] # make sdf function @@ -126,33 +206,55 @@ def sdf(invar, params, compute_sdf_derivatives=False): points = np.stack([invar["x"], invar["y"], invar["z"]], axis=1) # normalize triangles and points + store_triangles = np.array(triangles, dtype=np.float64) minx, maxx, miny, maxy, minz, maxz = _find_mins_maxs(points) max_dis = max(max((maxx - minx), (maxy - miny)), (maxz - minz)) - store_triangles = np.array(triangles, dtype=np.float64) - store_triangles[:, :, 0] -= minx - store_triangles[:, :, 1] -= miny - store_triangles[:, :, 2] -= minz - store_triangles *= 1 / max_dis - store_triangles = store_triangles.flatten() - points[:, 0] -= minx - points[:, 1] -= miny - points[:, 2] -= minz - points *= 1 / max_dis - points = points.astype(np.float64).flatten() + if self.pysdf_avalible: + minx, maxx, miny, maxy, minz, maxz = _find_mins_maxs(points) + max_dis = max(max((maxx - minx), (maxy - miny)), (maxz - minz)) + store_triangles = np.array(triangles, dtype=np.float64) + store_triangles[:, :, 0] -= minx + store_triangles[:, :, 1] -= miny + store_triangles[:, :, 2] -= minz + store_triangles *= 1 / max_dis + store_triangles = store_triangles.flatten() + + + points[:, 0] -= minx + points[:, 1] -= miny + points[:, 2] -= minz + points *= 1 / max_dis + points = points.astype(np.float64).flatten() + + # compute sdf values outputs = {} + mesh_new = copy(mesh) + mesh_new.vectors = store_triangles if airtight: - sdf_field, sdf_derivative = pysdf.signed_distance_field( - store_triangles, points, include_hit_points=True - ) - sdf_field = -np.expand_dims(max_dis * sdf_field, axis=1) + if self.pysdf_avalible: + sdf_field, sdf_derivative = pysdf.signed_distance_field( + store_triangles, points, include_hit_points=True + ) + else: + + # compute sdf values with mesh_to_sdf + vertices = mesh_new.vectors.reshape(-1, 3) + faces = np.arange(len(vertices)).reshape(-1, 3) + + trimesh_new = trimesh.Trimesh(vertices=vertices, + faces=faces) + sdf_field = mesh_to_sdf(trimesh_new, np.squeeze(points),surface_point_method='sample') + sdf_field = -np.expand_dims( sdf_field, axis=1) else: sdf_field = np.zeros_like(invar["x"]) outputs["sdf"] = sdf_field - + print('finished testelation') # get sdf derivatives if compute_sdf_derivatives: + if not self.pysdf_avalible: + raise Exception("sdf derivatives only curently work with pysdf") sdf_derivative = -(sdf_derivative - points) sdf_derivative = np.reshape( sdf_derivative, (sdf_derivative.shape[0] // 3, 3)