Mega Voxel Topology Optimization Example

Solving Massive Mega-Voxel Problems On The GPU

Here we demonstrate how you can use the multi-grid solver in FANTOM to solve mega-voxel TO problems in 3D. We will use the bridge problem again and use a structured mesh with 8M elements!

You can run this with 16GB of GPU VRAM since FANTOM’s multi-grid is matrix free at the highest level!

Importing The Necessary Tools

First we import the tools from pyFANTOM we will be using in this example.

from pyFANTOM.CUDA import(
    StructuredMesh3D, # Mesh Object for 3D Structured Meshes
    StructuredFilter3D, # Filter kernel for 3D Structured Meshes
    StructuredStiffnessKernel, # Stiffness Kernel for 3D Structured Meshes
    FiniteElement, # Finite Element Object to Setup Boundary Conditions and RHS
    CG, GMRES, MultiGrid, SPSOLVE, # Sparse Solvers
    MinimumCompliance, # TO Problem Definition
    PGD, MMA, OC # Nonlinear Optimizers
)

from pyFANTOM.Physics import LinearElasticity # Physics Model For FEA

# Other Imports
import numpy as np
import matplotlib.pyplot as plt
from tqdm import trange

Setting Up Physics and Geometry

The first step in using pyFANTOM is defining the domain of the problem and setting up the physics for FEA.

length = 1.0
width = 0.25
height = 0.25

# Create a physics model
physics = LinearElasticity(E=1.0, nu=1/3, thickness=1.0, type='PlaneStress')

# Create a mesh with 8M elements (512 x 128 x 128)
mesh = StructuredMesh3D(nx=512, ny=128, nz=128, lx=length, ly=height, lz=width, physics=physics)

Setting Up The Finite Element Solver

The next step in using pyFANTOM is setting up a stiffness kernel and finite-element boundary conditions. For mega-voxel problems, we use the MultiGrid solver with multiple levels to efficiently solve the large linear systems.

# Create a stiffness kernel
stiffness_kernel = StructuredStiffnessKernel(mesh = mesh)

# Create a solver
# For mega-voxel problems, MultiGrid with multiple levels is essential
solver = MultiGrid(
    mesh=mesh,
    kernel=stiffness_kernel,
    maxiter=100,
    tol=1e-4,
    n_level=6,  # 6 levels for the multi-grid hierarchy
    cycle='W',  # W-cycle for better convergence
    w_level=[2,3],  # Number of W-cycles at each level
    omega_boost=1.05  # Relaxation parameter boost
)

# Create a finite-element object
fe = FiniteElement(mesh=mesh, kernel=stiffness_kernel, solver=solver)

left_edge_nodes = np.where(mesh.nodes[:, 0] < 1e-12)[0]
right_edge_nodes = np.where(mesh.nodes[:, 0] > length - 1e-12)[0]

bottom_edge_nodes = np.where(mesh.nodes[:, 1] < 1e-12)[0]

fe.reset_dirichlet_boundary_conditions()
fe.reset_forces()
# Apply force to the bottom edge
fe.add_point_forces(
    node_ids=bottom_edge_nodes,
    positions=None,
    forces=np.array([[0.0, -1.0/len(bottom_edge_nodes), 0.0]])
)

# Apply zero displacement to the left edge and the right edge
fe.add_dirichlet_boundary_condition(
    node_ids=np.concatenate([left_edge_nodes, right_edge_nodes]), # Node IDs to apply the boundary condition to
    positions=None, # Alternative way to specify the boundary condition is to pass positions in space
    dofs=np.array([[1, 1, 1]]), # Degrees of freedom to apply the boundary condition to
    rhs=0.0 # Value to apply the boundary condition to
)

Setting Up TO Problem And Optimizer

Once the finite element is setup we can define the problem and optimizer to perform TO.

# Define the filter kernel for TO
filter_kernel = StructuredFilter3D(mesh=mesh, r_min=1.5)

# Define the TO problem
to_problem = MinimumCompliance(
    FE=fe,
    filter=filter_kernel,
    E_mul=[1.0], # You can pass a list of values to perform multi-material TO
    volume_fraction=[0.05], # You can pass a list of values to perform volume fraction control for each material
    void=1e-9, # You can pass a value to set void modulus
    penalty=3, # You can pass a value to set penalty factor
    # penalty_schedule = lambda p, i: (p-2)*np.round(3 * min(50, i) / 50)/3 + 2, # You can pass a function to set penalty schedule here
    heavyside= True, # You can pass a boolean to use heavyside or not
    eta=0.5, # You can pass a value to set eta - only used for heavyside
    beta=2.0, # You can pass a value to set beta, or a function of iteration to set beta schedule - only used for heavyside
)

# Define the optimizer
optimizer = PGD( # You can use any of the optimizers in pyFANTOM: OC, MMA, PGD
    problem=to_problem,
    change_tol=np.inf, # No change tolerance
    fun_tol=1e-4, # Function tolerance (convergence criterion for the optimizer)
)

Running The Optimization Loop

Now we can run the optimization loop and perform TO on the mega-voxel problem.

maximum_iterations = 300

Progressbar = trange(maximum_iterations, desc='Optimizer Iterations', leave=True)
objective_history = []
for i in Progressbar:
    optimizer.iter()

    Progressbar.set_postfix(
        optimizer.logs()
    )

    objective_history.append(optimizer.logs()['objective'])

    if optimizer.converged():
        print(f'Converged in {i} iterations')
        break

The mega-voxel optimization converged in 73 iterations with the following output:

Optimizer Iterations:  24%|██▍       | 73/300 [23:49<1:14:04, 19.58s/it, objective=244, variable change=5.01, function change=6.24e-6, iteration=75, residual=0.0135] Converged in 73 iterations

Now we can visualize the resulting topology:

to_problem.visualize_solution()

We can also verify the volume fraction:

(to_problem.desvars>0.5).sum()/to_problem.desvars.size

This returns approximately 0.05 (5%), confirming that the volume fraction constraint is satisfied.

Notes on Mega-Voxel Problems

• Memory Efficiency: FANTOM’s multi-grid solver is matrix-free at the highest level, allowing you to solve problems with 8M+ elements on GPUs with 16GB VRAM.

• Performance: The multi-grid solver with W-cycles provides excellent convergence rates even for very large problems, making it practical to solve mega-voxel topology optimization problems.

• Scalability: The solver scales well with the number of grid levels. For problems with 8M elements, using 6 levels provides a good balance between memory usage and convergence speed.