Smart Mesh

notebooks/006_dipolarWakefield_mesh_refinement.ipynb

@@ -116,7 +116,8 @@
     "grid = GridFIT3D(xmin, xmax, ymin, ymax, zmin, zmax, Nx, Ny, Nz,\n",
     "                stl_solids=stl_solids,\n",
     "                stl_materials=stl_materials,\n",
-    "                use_mesh_refinement=False,\n",
+    "                use_mesh_refinement=True,\n",
+    "                snap_tol=1e-4,\n",
     "                refinement_tol=1e-8,)"
    ]
   },
@@ -130,6 +131,16 @@
     "# grid.inspect()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# [OPTIONAL] Look at Snap Points\n",
+    "# grid.plot_snap_points(snap_tol=1e-4)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,

tests/test_005_pml_boundaries.py

@@ -48,7 +48,7 @@ def test_reflection_planewave(self):
         solver.dt = 1/f/200 #ensure right amplitude
 
         # Simulation
-        Nt = int(2.0*(zmax-zmin-solver.n_pml*solver.dz)/c/solver.dt)
+        Nt = int(2.0*(zmax-zmin-solver.n_pml*np.min(solver.dz))/c/solver.dt) #only temporary solution, it takes the smallest dz as measure for pml time
 
         for n in tqdm(range(Nt)):
             planeWave.update(solver, n*solver.dt)

tests/test_007_mpi_lossy_cavity.py

@@ -320,8 +320,8 @@ def test_long_impedance(self):
             assert np.allclose(np.real(wake.Z)[::20], np.real(self.Z), rtol=0.1), "Real Impedance samples failed"
             assert np.allclose(np.imag(wake.Z)[::20], np.imag(self.Z), rtol=0.1), "Imag Impedance samples failed"
             assert np.cumsum(np.abs(wake.Z))[-1] == pytest.approx(250910.51090497518, 0.1), "Abs Impedance cumsum failed"
-
-    def test_log_file(self):
+    
+    def test_log_file(self):       
         # Helper function to compare nested dicts with float tolerance
         def assert_dict_allclose(d1, d2, rtol=1e-6, atol=1e-12, path=""):
             assert set(d1.keys()) == set(d2.keys()), \
@@ -337,11 +337,14 @@ def assert_dict_allclose(d1, d2, rtol=1e-6, atol=1e-12, path=""):
 
                 # floats
                 elif isinstance(v1, float) and isinstance(v2, float):
-                    assert np.isclose(v1, v2, rtol=rtol, atol=atol), \
-                        f"Float mismatch at {p}: {v1} != {v2}"
+                    if k == 'dt':
+                        assert v1 <=v2, f"Timestep bigger than for uniform grid"
+                    else:
+                        assert np.isclose(v1, v2, rtol=rtol, atol=atol), \
+                            f"Float mismatch at {p}: {v1} != {v2}"
 
                 # lists/tuples/arrays
-                elif isinstance(v1, (list, tuple)) and isinstance(v2, (list, tuple)):
+                elif isinstance(v1, (list, tuple, np.ndarray)) and isinstance(v2, (list, tuple, np.ndarray)):
                     assert len(v1) == len(v2), f"Length mismatch at {p}"
                     for i, (a, b) in enumerate(zip(v1, v2)):
                         if isinstance(a, float) and isinstance(b, float):
@@ -350,6 +353,17 @@ def assert_dict_allclose(d1, d2, rtol=1e-6, atol=1e-12, path=""):
                         else:
                             assert a == b, f"Value mismatch at {p}[{i}]: {a} != {b}"
 
+                # list vs single float
+                elif isinstance(v1, (list, tuple, np.ndarray)) and isinstance(v2, float):
+                    for i, a in enumerate(v1):
+                        assert np.isclose(a, v2, rtol=rtol, atol=atol), \
+                            f"Float mismatch at {p}[{i}]: {a} != {v2}"
+
+                elif isinstance(v1, float) and isinstance(v2, (list, tuple, np.ndarray)):
+                    for i, b in enumerate(v2):
+                        assert np.isclose(v1, b, rtol=rtol, atol=atol), \
+                            f"Float mismatch at {p}[{i}]: {v1} != {b}"
+
                 # everything else → exact match
                 else:
                     assert v1 == v2, f"Mismatch at {p}: {v1} != {v2}"
@@ -368,4 +382,4 @@ def assert_dict_allclose(d1, d2, rtol=1e-6, atol=1e-12, path=""):
         # Compare log dict contents
         assert_dict_allclose(solver.logger.grid, self.gridLogs)
         assert_dict_allclose(solver.logger.solver, self.solverLogs)
-        assert_dict_allclose(solver.logger.wakeSolver, self.wakeSolverLogs)
+        assert_dict_allclose(solver.logger.wakeSolver, self.wakeSolverLogs)

tests/test_011_smart_mesh.py

@@ -0,0 +1,110 @@
+import os, sys
+import numpy as np
+import pyvista as pv
+
+sys.path.append('../wakis')
+
+from wakis import SolverFIT3D
+from wakis import GridFIT3D 
+from wakis import WakeSolver
+
+import pytest 
+
+@pytest.mark.slow
+class TestSmartMesh:
+
+    zedges = [-0.05,-0.0075,0.0075,0.05]
+    xedges = [-0.025,-0.005, 0.005, 0.025]
+    yedges = [-0.025,-0.005, 0.005, 0.025]
+
+    print("\n---------- Initializing simulation ------------------")
+    # Number of mesh cells
+    Nx = 50
+    Ny = 50
+    Nz = 80
+
+    # Embedded boundaries
+    stl_file = 'tests/stl/001_cubic_cavity.stl' 
+    surf = pv.read(stl_file)
+
+    stl_solids = {'cavity': stl_file}
+    stl_materials = {'cavity': 'vacuum'}
+
+    # Domain bounds
+    xmin, xmax, ymin, ymax, zmin, zmax = surf.bounds
+
+    refinement_tol=1e-8
+    snap_tol=1e-2
+    # set grid and geometry
+    global grid
+    grid = GridFIT3D(xmin, xmax, ymin, ymax, zmin, zmax, Nx, Ny, Nz, 
+                    stl_solids=stl_solids, 
+                    stl_materials=stl_materials,
+                    use_mesh_refinement=True,
+                    snap_tol=snap_tol,
+                    refinement_tol=refinement_tol)
+
+    # Beam parameters
+    beta = 1.          # beam beta
+    sigmaz = 18.5e-3*beta    #[m]
+    q = 1e-9            #[C]
+    xs = 0.             # x source position [m]
+    ys = 0.             # y source position [m]
+    xt = 0.             # x test position [m]
+    yt = 0.             # y test position [m]
+
+    global wake
+    skip_cells = 12  # no. cells to skip in WP integration
+    wake = WakeSolver(q=q, sigmaz=sigmaz, beta=beta,
+                xsource=xs, ysource=ys, xtest=xt, ytest=yt,
+                save=False, Ez_file='tests/011_Ez.h5',
+                skip_cells=skip_cells,
+                )
+
+    # boundary conditions
+    bc_low=['pec', 'pec', 'pec']
+    bc_high=['pec', 'pec', 'pec']
+
+    # set Solver object
+    solver = SolverFIT3D(grid, wake,
+                        bc_low=bc_low, bc_high=bc_high,
+                        use_stl=True, bg='pec')
+
+    wakelength = 1. #[m]
+    solver.wakesolve(wakelength=wakelength, save_J=False)
+    os.remove('tests/011_Ez.h5')
+
+    def test_grid_generation(self):
+        global grid
+        for edg in self.zedges:
+            diff = np.min(np.abs(grid.z-edg))
+            assert diff <= 1e-8 + 1e-2, "Mesh not inside the tolerance at the z edges"
+
+        for edg in self.yedges:
+            diff = np.min(np.abs(grid.y-edg))
+            assert diff <= 1e-8 + 1e-2, "Mesh not inside the tolerance at the y edges"
+
+        for edg in self.xedges:
+            diff = np.min(np.abs(grid.x-edg))
+            assert diff <= 1e-8 + 1e-2, "Mesh not inside the tolerance at the x edges"
+
+        zdiff=np.abs(np.diff(grid.dz))
+        for cell in range(len(grid.dz)-1):
+            assert zdiff[cell]<=0.2*grid.dz[cell], "Mesh variance in z is too big"
+
+        ydiff=np.abs(np.diff(grid.dy))
+        for cell in range(len(grid.dy)-1):
+            assert ydiff[cell]<=0.2*grid.dy[cell], "Mesh variance in y is too big"
+
+        xdiff=np.abs(np.diff(grid.dx))
+        for cell in range(len(grid.dx)-1):
+            assert xdiff[cell]<=0.2*grid.dx[cell], "Mesh variance in x is too big"
+
+        assert np.min(grid.dz) > 0.25 * (grid.zmax-grid.zmin)/grid.Nz, "Smallest z difference is too small"
+        assert np.min(grid.dy) > 0.25 * (grid.ymax-grid.ymin)/grid.Ny, "Smallest y difference is too small"
+        assert np.min(grid.dx) > 0.25 * (grid.xmax-grid.xmin)/grid.Nx, "Smallest x difference is too small"
+
+    def test_nonuniform_simulation(self):
+        global wake
+        assert np.abs(wake.Z[0]) < 0.01 * np.abs(np.max(wake.Z)), "Charge Accumulation - DC component"
+        # TODO: Test that checks the continuity equation

wakis/gridFIT3D.py

@@ -9,6 +9,7 @@
 import h5py
 
 from scipy.optimize import least_squares
+import time
 
 from .field import Field
 from .logger import Logger
@@ -129,10 +130,9 @@ def __init__(self, xmin=None, xmax=None,
             )
 
         #TODO: allow non uniform dx, dy, dz
-        self.dx = np.min(np.diff(self.x))
-        self.dy = np.min(np.diff(self.y))
-        #self.dz = np.min(np.diff(self.z))
-        self.dz = (self.zmax - self.zmin)/self.Nz
+        self.dx = np.diff(self.x)
+        self.dy = np.diff(self.y)
+        self.dz = np.diff(self.z)
         self.update_logger(['Nx', 'Ny', 'Nz', 'dx', 'dy', 'dz'])
         self.update_logger(['xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax'])
 
@@ -216,9 +216,12 @@ def _compute_grid(self):
         self.L.field_z = Z[1:, 1:, 1:] - Z[:-1, :-1, :-1]
 
         self.iA = Field(self.Nx, self.Ny, self.Nz)
-        self.iA.field_x = np.divide(1.0, self.L.field_y * self.L.field_z)
-        self.iA.field_y = np.divide(1.0, self.L.field_x * self.L.field_z)
-        self.iA.field_z = np.divide(1.0, self.L.field_x * self.L.field_y)
+        aux = self.L.field_y * self.L.field_z
+        self.iA.field_x = np.divide(1.0, aux , out=np.zeros_like(aux), where=aux!=0)
+        aux = self.L.field_x * self.L.field_z
+        self.iA.field_y = np.divide(1.0, aux , out=np.zeros_like(aux), where=aux!=0)
+        aux = self.L.field_x * self.L.field_y
+        self.iA.field_z = np.divide(1.0, aux , out=np.zeros_like(aux), where=aux!=0)
 
         # tilde grid ~G
         self.tx = (self.x[1:]+self.x[:-1])/2
@@ -288,6 +291,7 @@ def _mpi_initialize(self):
             self.Nz += self.n_ghosts
 
         self.z = np.linspace(self.zmin, self.zmax, self.Nz+1)
+        self.dz = np.diff(self.z) # only uniform grid possible with MPI
 
     def mpi_gather_asGrid(self):
         _grid = None
@@ -357,7 +361,7 @@ def _prepare_stl_dicts(self):
 
     def _mark_cells_in_stl(self):
         # Obtain masks with grid cells inside each stl solid
-        stl_tolerance = np.min([self.dx, self.dy, self.dz])*self.stl_tol
+        stl_tolerance = np.min([np.min(self.dx), np.min(self.dy), np.min(self.dz)])*self.stl_tol
         progress_bar = False
         if self.Nx*self.Ny*self.Nz > 5e6 and self.verbose:
             progress_bar = True
@@ -418,6 +422,10 @@ def _compute_snap_points(self, snap_solids=None, snap_tol=1e-8):
                 model = model + solid
 
         edges = model.extract_feature_edges(boundary_edges=True, manifold_edges=False)
+        print(edges)
+        print("min abs(Y):", np.min(np.abs(edges.points[:,1])))
+        print("min abs(X):", np.min(np.abs(edges.points[:,0])))
+        print("min abs(Z):", np.min(np.abs(edges.points[:,2])))
 
         # Extract points lying in the X-Z plane (Y ≈ 0)
         xz_plane_points = edges.points[np.abs(edges.points[:, 1]) < snap_tol]
@@ -566,12 +574,11 @@ def _refine_xyz_axis(self, method='insert', tol=1e-6):
         self.Nx = len(self.x) - 1
         self.Ny = len(self.y) - 1
         self.Nz = len(self.z) - 1
-        self.dx = np.min(np.diff(self.x))  #TODO: should this be an array?
-        self.dy = np.min(np.diff(self.y))
-        self.dz = np.min(np.diff(self.z))
+        self.dx = np.diff(self.x)
+        self.dy = np.diff(self.y)
+        self.dz = np.diff(self.z)
 
-        if self.verbose:
-            print(f"Refined grid: Nx = {len(self.x)}, Ny ={len(self.y)}, Nz = {len(self.z)}")
+        print(f"Refined grid: Nx = {self.Nx}, Ny ={self.Ny}, Nz = {self.Nz}")
 
     def _assign_colors(self):
         '''Classify colors assigned to each solid
@@ -1016,9 +1023,9 @@ def load_from_h5(self, filename):
         self.Nx = len(self.x) - 1
         self.Ny = len(self.y) - 1
         self.Nz = len(self.z) - 1
-        self.dx = np.min(np.diff(self.x))
-        self.dy = np.min(np.diff(self.y))
-        self.dz = np.min(np.diff(self.z))
+        self.dx = np.diff(self.x)
+        self.dy = np.diff(self.y)
+        self.dz = np.diff(self.z)
         self.xmin, self.xmax = self.x[0], self.x[-1]
         self.ymin, self.ymax = self.y[0], self.y[-1]
         self.zmin, self.zmax = self.z[0], self.z[-1]

wakis/routines.py

@@ -85,6 +85,7 @@ def emsolve(self, Nt, source=None, callback=None,
 
             for hf in hfs:
                 hf['x'], hf['y'], hf['z'] = self.x, self.y, self.z
+                hf['dx'], hf['dy'], hf['dz'] = self.grid.dx, self.grid.dy, self.grid.dz
                 hf['t'] = np.arange(0, Nt*self.dt, save_every*self.dt)
 
             if subdomain is not None:
@@ -247,6 +248,7 @@ def wakesolve(self, wakelength,
         self.v = self.wake.v
         if self.use_mpi:  #E- should it be zmin, zmax instead?
             z = self.Z # use global coords
+            dz = np.diff(self.Z)
             zz = slice(0, self.NZ)
         else:
             z = self.z
@@ -270,20 +272,24 @@ def wakesolve(self, wakelength,
             if self.rank==0:
                 hf = h5py.File(self.Ez_file, 'w')
                 hf['x'], hf['y'], hf['z'] = self.x[xx], self.y[yy], z[zz]
+                hf['dx'], hf['dy'], hf['dz'] = self.grid.dx, self.grid.dy, dz
                 hf['t'] = np.arange(0, Nt*self.dt, self.dt)
 
                 if save_J:
                     hfJ = h5py.File('Jz.h5', 'w')
                     hfJ['x'], hfJ['y'], hfJ['z'] = self.x[xx], self.y[yy], z[zz]
+                    hfJ['dx'], hfJ['dy'], hfJ['dz'] = self.grid.dx, self.grid.dy, dz
                     hfJ['t'] = np.arange(0, Nt*self.dt, self.dt)
         else:
             hf = h5py.File(self.Ez_file, 'w')
             hf['x'], hf['y'], hf['z'] = self.x[xx], self.y[yy], z[zz]
+            hf['dx'], hf['dy'], hf['dz'] = self.grid.dx, self.grid.dy, self.grid.dz
             hf['t'] = np.arange(0, Nt*self.dt, self.dt)
 
             if save_J:
                 hfJ = h5py.File('Jz.h5', 'w')
                 hfJ['x'], hfJ['y'], hfJ['z'] = self.x[xx], self.y[yy], z[zz]
+                hfJ['dx'], hfJ['dy'], hfJ['dz'] = self.grid.dx, self.grid.dy, self.grid.dz
                 hfJ['t'] = np.arange(0, Nt*self.dt, self.dt)
 
         def save_to_h5(self, hf, field, x, y, z, comp, n):

wakis/solverFIT3D.py

@@ -176,6 +176,9 @@ def __init__(self, grid, wake=None, cfln=0.5, dt=None,
                             hstack([self.Pz, sparse_mat((N,N)), -self.Px]),
                             hstack([-self.Py, self.Px, sparse_mat((N,N))])
                         ], dtype=np.int8)
+
+        self.S = hstack([self.Px, self.Py, self.Pz])
+        print(self.S*self.C)
 
         # Boundaries
         if verbose:
@@ -220,7 +223,7 @@ def __init__(self, grid, wake=None, cfln=0.5, dt=None,
             print('Calculating maximal stable timestep...')
         self.cfln = cfln
         if dt is None:
-            self.dt = cfln / (c_light * np.sqrt(1 / self.grid.dx ** 2 + 1 / self.grid.dy ** 2 + 1 / self.grid.dz ** 2))
+            self.dt = cfln / (c_light * np.sqrt(1 / np.min(self.grid.dx) ** 2 + 1 / np.min(self.grid.dy) ** 2 + 1 / np.min(self.grid.dz) ** 2))
         else:
             self.dt = dt
         self.dt = dtype(self.dt)
@@ -244,8 +247,8 @@ def __init__(self, grid, wake=None, cfln=0.5, dt=None,
         self.iDmu = diags(self.imu.toarray(), shape=(3*N, 3*N), dtype=self.dtype)
         self.Dsigma = diags(self.sigma.toarray(), shape=(3*N, 3*N), dtype=self.dtype)
 
-        self.tDsiDmuiDaC = self.tDs * self.iDmu * self.iDa * self.C
-        self.itDaiDepsDstC = self.itDa * self.iDeps * self.Ds * self.C.transpose()
+        self.tDsiDmuiDaC = self.iDa * self.iDmu * self.C * self.Ds
+        self.itDaiDepsDstC = self.iDeps * self.itDa * self.C.transpose() * self.tDs
 
         if imported_mkl and not self.use_gpu: # MKL backend for CPU
             if verbose:
@@ -272,6 +275,9 @@ def __init__(self, grid, wake=None, cfln=0.5, dt=None,
         self.solverInitializationTime = time.time() - t0
         self.update_logger(['solverInitializationTime'])
 
+        self.solverInitializationTime = time.time() - t0
+        self.update_logger(['solverInitializationTime'])
+
     def update_tensors(self, tensor='all'):
         '''Update tensor matrices after
         Field ieps, imu or sigma have been modified
@@ -300,8 +306,8 @@ def update_tensors(self, tensor='all'):
 
         if self.verbose:
             print('Re-Pre-computing ...')
-        self.tDsiDmuiDaC = self.tDs * self.iDmu * self.iDa * self.C
-        self.itDaiDepsDstC = self.itDa * self.iDeps * self.Ds * self.C.transpose()
+        self.tDsiDmuiDaC = self.iDa * self.iDmu * self.C * self.Ds
+        self.itDaiDepsDstC = self.iDeps * self.itDa * self.C.transpose() * self.tDs
         self.step_0 = False
 
     def _one_step(self):
@@ -322,7 +328,8 @@ def _one_step(self):
 
         #include current computation
         if self.use_conductivity:
-            self.J.fromarray(self.sigma.toarray()*self.E.toarray())
+            self.J.fromarray(self.sigma.toarray()*self.E.toarray()
+                             )
 
     def one_step_mkl(self):
         if self.step_0: