We present and distribute a FreeFem++ Toolbox for the parallel computing of two- or three-dimensional liquid–solid phase-change systems involving natural convection. FreeFem++ (www.freefem.org) is a free finite-element software available for all existing operating systems. We use the recent library ffddm that makes available in FreeFem++ state-of-the-art scalable Schwarz domain decomposition methods (DDM). The single domain approach used in our previous contribution (Rakotondrandisa et al., 2020) is adapted for the use of the DDM method. As a result, the computational time is considerably reduced for 2D configurations and furthermore 3D problems become affordable. The numerical method is based on an enthalpy-porosity model. The same set of equations is solved in both liquid and solid phases: the incompressible Navier–Stokes equations with Boussinesq approximation for thermal effects. A Carman–Kozeny-type penalty term is added to the momentum equations to bring progressively the velocity to zero into the solid. Model equations are discretized using Galerkin triangular or tetrahedral finite elements. The coupled system of equations is integrated in time using a second-order Gear implicit scheme. The resulting discrete equations are solved using a Newton algorithm. The DDM approach is based on an overlapping Schwarz method. The mesh is first split in subdomains using Scotch or Metis libraries. The final linear system is then solved in parallel using a GMRES Krylov method, with a Restricted Additive Schwarz (RAS) preconditioner. The mesh is adapted during the computation using metrics control. The 3D-mesh adaptivity uses the mmg (www.mmgtools.org) open source library. Parallel 2D and 3D computations of benchmark cases of increasing difficulty are presented: natural convection of air, natural convection of water, melting or solidification of a phase-change material, and, finally, a water freezing case. For each case, careful validations are provided and the performance of the code is assessed. The robustness of the Toolbox in 3D is also demonstrated by adapting the number of processors to the number of tetrahedra, which can considerably vary after the mesh adaptation. Program summaryProgram Title: PCM_Toolbox_DDM_2D and PCM_Toolbox_DDM_3DCPC Library link to program files:http://dx.doi.org/10.17632/dk49rfrz9y.1Licensing provisions: Apache License, 2.0Programming language: FreeFem++(www.freefem.org), mmg (www.mmgtools.org)Nature of problem: The software is scoped to parallel computations of 2D or 3D configurations of liquid–solid phase-change problems with convection in the liquid phase. Natural convection, melting and solidification processes are illustrated in the paper. The software can be easily modified to take into account different related physical models.Solution method: We use a single domain approach, solving the incompressible Navier–Stokes equations with Boussinesq approximation in both liquid and solid phases. A Carman–Kozeny-type penalty term is added to the momentum equations to bring the velocity to zero into the solid phase. An enthalpy model is used in the energy equation to take into account the phase change. Discontinuous variables (latent heat, material properties) are regularized through an intermediate (mushy) region. Space discretization is based on Galerkin triangular/tetrahedral finite elements. A second order Gear implicit scheme is used for the time integration of the coupled system of equations. The resulting discrete equations are solved using a Newton algorithm. Piecewise quadratic (P2) finite-elements are used for the velocity and piecewise linear (P1) for the pressure. For the temperature both P2 or P1 discretization are possible. The mesh is first split in subdomains using Scotch or Metis libraries, which are interfaced with FreeFem++. Then, a Schwarz domain decomposition method is used through the FreeFem++library ffddm. The final linear system is solved in parallel using a GMRES Krylov method, with a Restricted Additive Schwarz (RAS) preconditioner. Mesh adaptivity using metrics control makes possible the optimization of the distribution of mesh elements. For 3D case, the mmg open source library is used to adapt the mesh.