Porous materials with hierarchical structures are widely used in various engineering and scientific applications due to their enhanced properties. However, the hierarchical nature of their structures presents challenges when it comes to modeling. In this work, we propose a multilevel numerical homogenization method based on the concept of splitting the Representative Volume Element (RVE) into subdomains. A random porous structure is modeled at the sub-micron and micron levels to determine its apparent modulus tensor using the previously presented framework developed in FEniCSx. To provide better separation of scale levels and capture variability in the structure, higher-order homogenization is implemented at the sub-micron and micron scales. After this, using the obtained properties' distributions, a dense heterogeneous media was modeled with further performing first-order homogenization to determine the effective tensor of elastic moduli at the mesoscale level. We test the proposed framework for various boundary value problems to establish the effective modulus tensor and compare the results with standard first-order homogenization for a full-scale model of the same problem. The results are in good agreement with the classical method but at significantly lower computational costs, which is the promise of the proposed method.