Note

Go to the end to download the full example code.

Multigrid preconditioned CG solver

Significant speedups can be achieved using CG with a multigrid preconditioner to solve the linear system of equations

The speedup enables large-scale topology optimization problems, which is achieved by projecting the finite element problem onto a series of coarser grids, and only factorizing the finite element matrix at the coarsest grid. The solution is then interpolated back to the fine grid. This process is used as preconditioner for the CG iterations, enabling fast convergence and optimization in 3D.

This example is based on Compliance minimization and only differs in the solver used. The user has full control over the number of multigrid levels and their internal settings (see pymoto.solvers.GeometricMultigrid and pymoto.solvers.CG for more details).

References:

Amir, O., Aage, N., & Lazarov, B. S. (2014). On multigrid-CG for efficient topology optimization. Structural and Multidisciplinary Optimization, 49(5), 815-829. DOI: https://doi.org/10.1007/s00158-013-1015-5

 21 import numpy as np
 22 import pymoto as pym
 23
 24 # nx, ny, nz = 32, 64, 64  # Domain sizes (# elements). Set nz to zero for the 2D problem, nz > 0 runs a 3D problem
 25 nx, ny, nz = 512, 256, 0  # 2D
 26 xmin = 1e-9
 27 filter_radius = 2.0
 28 volfrac = 0.5
 29 thermal = False  # True = Static thermal analysis; False = Static mechanical analysis will be done
 30
 31 if __name__ == "__main__":
 32     print(__doc__)
 33
 34     if nz == 0:  # 2D analysis
 35         # Generate a grid
 36         domain = pym.VoxelDomain(nx, ny)
 37
 38         if thermal:
 39             ndof = 1  # Number of dofs per node
 40             # Get dof numbers at the boundary
 41             boundary_dofs = domain.nodes[0, ny//4:((ny + 1) - ny // 4)].flatten()
 42
 43             # Make a force vector
 44             force_dofs = domain.nodes[1:, :].flatten()
 45
 46         else:  # Mechanical
 47             ndof = 2
 48             # Calculate boundary dof indices
 49             boundary_nodes = domain.nodes[0, :].flatten()
 50             boundary_dofs = domain.get_dofnumber(boundary_nodes, ndof=ndof)
 51
 52             # Which dofs to put a force on? The 1 is added for a force in y-direction (x-direction would be zero)
 53             force_dofs = ndof * domain.nodes[nx, ny//2] + 1
 54
 55     else:
 56         domain = pym.VoxelDomain(nx, ny, nz)
 57         boundary_nodes = domain.nodes[0, :, :].flatten()
 58
 59         if thermal:
 60             boundary_dofs = boundary_nodes
 61             force_dofs = domain.nodes[1:, :, :].flatten()
 62             ndof = 1
 63
 64         else:  # Mechanical
 65             ndof = 3
 66             boundary_dofs = domain.get_dofnumber(boundary_nodes, ndof=ndof)
 67             force_dofs = ndof * domain.get_nodenumber(nx, ny // 2, nz // 2) + 2  # Z-direction
 68
 69     if domain.nnodes > 1e+6:
 70         print("Too many nodes :(")  # Safety, to prevent overloading the memory in your machine
 71         exit()
 72
 73     # Generate a force vector
 74     f = np.zeros(domain.nnodes * ndof)
 75     f[force_dofs] = 1.0  # Uniform force of 1.0 at all selected dofs
 76
 77     # Make design vector and fill with initial values
 78     sx = pym.Signal('x', state=np.ones(domain.nel) * volfrac)
 79
 80     # Start building the modular network
 81     with pym.Network(print_timing=True) as func:
 82         # Filter
 83         sxfilt = pym.FilterConv(domain, radius=filter_radius)(sx)
 84         sx_analysis = sxfilt
 85         sx_analysis.tag = 'design'
 86
 87         # Show the design on the screen as it optimizes (plotting domain in 3D may take significant time)
 88         if nz == 0:
 89             pym.PlotDomain(domain, saveto="out/design", clim=[0, 1])(sx_analysis)
 90
 91         # SIMP material interpolation
 92         sSIMP = pym.MathExpression(f"{xmin} + {1.0 - xmin}*inp0^3")(sx_analysis)
 93
 94         # System matrix assembly module
 95         if thermal:
 96             sK = pym.AssemblePoisson(domain, bc=boundary_dofs)(sSIMP)
 97         else:
 98             sK = pym.AssembleStiffness(domain, bc=boundary_dofs)(sSIMP)
 99
100         # Linear system solver. The linear solver can be chosen by uncommenting any of the following lines.
101         solver = None  # Default (solver is automatically chosen based on matrix properties)
102         # solver = pym.solvers.SolverSparsePardiso()  # Requires Intel MKL installed
103         # solver = pym.solvers.SolverSparseCholeskyCVXOPT()  # Requires cvxopt installed
104         # solver = pym.solvers.SolverSparseCholeskyScikit()  # Requires scikit installed
105
106         ''' Iterative solver: CG with geometric multigrid preconditioning '''
107         # Set up fine level multigrid operator
108         mg1 = pym.solvers.GeometricMultigrid(domain)
109         mgs = [mg1] # List of multigrid levels, starting with the finest grid
110         while True:  # Setup multiple sub-levels
111             # Stop if any of the new dimensions is odd
112             if any([nn % 2 != 0 for nn in mgs[-1].sub_domain.size]):
113                 break
114
115             # Stop if any of the coarse dimensions is smaller than 8
116             if any(mgs[-1].sub_domain.size < 8):
117                 break
118
119             # Create a new coarser grid
120             mgs.append(pym.solvers.GeometricMultigrid(mgs[-1].sub_domain))
121             mgs[-2].inner_level = mgs[-1]  # Set the inner level to the next coarser grid
122
123         print(f"Multigrid levels: {len(mgs)}")
124         print(f"Coarsest grid size: {mgs[-1].sub_domain.nelx} x {mgs[-1].sub_domain.nely} " +
125               (f"x {mgs[-1].sub_domain.nelz}" if nz > 0 else ""))
126
127         # Set up the solver (comment out to use the default factorization, try this to see the difference in time)
128         solver = pym.solvers.CG(preconditioner=mg1, verbosity=1, tol=1e-5)
129
130         ''' From here the rest of the example is identical to ex_compliance.py '''
131         su = pym.LinSolve(hermitian=True, solver=solver)(sK, f)
132
133         # Output the design, deformation, and force field to a Paraview file
134         pym.WriteToVTI(domain, saveto='out/dat.vti')(sx_analysis, su, f)
135
136         # Compliance calculation c = f^T u
137         scompl = pym.EinSum('i,i->')(su, f)
138         scompl.tag = 'compliance'
139
140         # MMA needs correct scaling of the objective
141         sg0 = pym.Scaling(scaling=100.0)(scompl)
142         sg0.tag = "objective"
143
144         # Calculate the volume of the domain by adding all design densities together
145         svol = pym.EinSum('i->')(sx_analysis)
146         svol.tag = 'volume'
147
148         # Volume constraint
149         sg1 = pym.MathExpression(f'10*(inp0/{domain.nel} - {volfrac})')(svol)
150         sg1.tag = "volume constraint"
151
152     # Maybe you want to check the design-sensitivities?
153     do_finite_difference = False
154     if do_finite_difference:
155         pym.finite_difference(sx, [sg0, sg1], function=func, dx=1e-4)
156         exit()
157
158     with func:
159         pym.PlotIter()(sg0, sg1)  # Plot iteration history
160
161     # Do the optimization with MMA
162     pym.minimize_mma([sx], [sg0, sg1], function=func)
163
164     # Do the optimization with OC
165     # pym.minimize_oc(sx, sg0, function=func)
166
167     # Here you can do some post processing
168     print("The optimization has finished!")
169     print(f"The final compliance value obtained is {scompl.state}")
170     print(f"The maximum {'temperature' if thermal else 'displacement'} is {max(np.absolute(su.state))}")

Gallery generated by Sphinx-Gallery