Selected Publications

ACM Transactions on Graphics, v. 44, n. 4, 2025.

Filipe Nascimento, Fabricio S. Sousa, Afonso Paiva  

Abstract: 

Powder-snow avalanches are natural phenomena that result from an instability in the snow cover on a mountain relief. It begins with a dense avalanche core moving fast down the mountain. During its evolution, the snow particles in the avalanche front mix with the air, forming a suspended turbulent cloud of snow dust surrounding the dense snow avalanche. This paper introduces a physically-based framework using the Finite Volume Method to simulate powder-snow avalanches under complex terrains. Specifically, the primary goal is to simulate the turbulent snow cloud dynamics within the avalanche in a visually realistic manner. Our approach relies on a multi-layer model that splits the avalanche into two main layers: dense and powder-snow. The dense-snow layer flow is simulated by solving a type of Shallow Water Equations suited for intricate basal surfaces, known as the Savage-Hutter model. The powder-snow layer flow is modeled as a two-phase mixture of miscible fluids and simulated using Navier-Stokes equations. Moreover, we propose a novel model for the transition layer, which is responsible for coupling the avalanche main layers, including the snow mass injected into the powder-snow cloud from the snow entrainment processes and its injection velocity. In brief, our framework comprehensively simulates powder-snow avalanches, allowing us to render convincing animations of one of the most complex gravity-driven flows.

DOI Open Access

Ocean Engineering, v. 296, 116699, 2024.

Lucas S. Pereira, Rubens A. Amaro Jr., Liang-Yee Cheng, Fabricio S. Sousa, Gustavo M. Karuka  

Abstract: 

MPS-VG, a Virtual Grating (VG) model for the Lagrangian mesh-free Moving Particle Semi-implicit (MPS) method is proposed for replacing conventional particle-based solid modeling of gratings with a set of thin inclined slats. Unlike most approaches for perforated wave energy dampening devices, in which the flow through the device is simplified by pressure loss or damping effects without flow deflection, MPS-VG models the angular deviation caused by hydrodynamic impact on inclined slats. Both accuracy and computational performance of the model were checked through a simulation of wet dam break scenarios with the grating structures placed horizontally or vertically. The results were compared with those from fully particle-based modeling. MPS-VG correctly predicted complex wave-structure interactions using a relatively low-resolution model and significantly reduced processing time and memory storage compared to conventional particle-based MPS modeling. The evaluation of the performance of the gratings with inclined slats as wave energy dampers revealed the horizontal gratings outperformed the vertical ones. Therefore, qualitative and quantitative agreements strengthened the potential of MPS-VG as a practical and computationally efficient tool for the study of multi-scale phenomena of wave impacts on grating with inclined slats.

DOI

Journal of Computational Physics, v. 473, p. 111681, 2023.

E. Abreu, P. Ferraz, A. M. Espírito Santo, F. Pereira, L. G. C. Santos, Fabricio S. Sousa  

Abstract: 

Multiscale methods for second order elliptic equations based on non-overlapping domain decomposition schemes have great potential to take advantage of multi-core, state-of-the-art parallel computers. These methods typically involve solving local boundary value problems followed by the solution of a global interface problem. Known iterative procedures for the solution of the interface problem have typically slow convergence, increasing the overall cost of the multiscale solver. To overcome this problem we develop a scalable recursive solution method for such interface problem that replaces the global problem by a family of small interface systems associated with adjacent subdomains, in a hierarchy of nested subdomains. Then, we propose a novel parallel algorithm to implement our recursive formulation in multi-core devices using the Multiscale Robin Coupled Method by Guiraldello et al. (2018), that can be seen as a generalization of several multiscale mixed methods. Through several numerical studies we show that the new algorithm is very fast and exhibits excellent strong and weak scalability. We consider very large problems, that can have billions of discretization cells, motivated by the numerical simulation of subsurface flows.

DOI | Pre-print 

Journal of Computational Science, v. 60, p. 101592, 2022.

Franciane F. Rocha, Fabricio S. Sousa, Roberto F. Ausas, Gustavo C. Buscaglia, Felipe Pereira

Abstract: 

In the presence of strong heterogeneities, it is well known that the use of explicit schemes for the transport of species in a porous medium suffers from severe restrictions on the time step. This has led to the development of implicit schemes that are increasingly favored by practitioners for their computational efficiency. The transport equation requires knowledge of the velocity field, which results from an elliptic problem (Darcy problem) that is the most expensive part of the computation. When considering large reservoirs, a cost-effective way of approximating the Darcy problems is using multiscale domain decomposition (MDD) methods. They allow for the pressure and velocity fields to be computed on coarse meshes (large scale), while detailed basis functions are defined locally, usually in parallel, in a much finer grid (small scale). In this work we adopt the Multiscale Robin Coupled Method (MRCM, Guiraldello et al., 2018; Rocha et al., 2020), which is a generalization of previous MDD methods that allows for great flexibility in the choice of interface spaces. In this article we investigate the combination of the MRCM with implicit transport schemes. A sequentially implicit strategy is proposed, with different trust-region algorithms ensuring the convergence of the transport solver. The method is assessed on several very stringent 2D two-phase problems, demonstrating its stability even for large time steps. It is also shown that the best accuracy is achieved by considering recently introduced non-polynomial interface spaces, since polynomial spaces are not optimal for high-contrast channelized permeability fields.

DOI Open Access | PDF 

Computational Geosciences, v. 26, p. 481-501, 2022.

Alfredo Jaramillo, Rafael T. Guiraldello, Stevens Paz, Roberto F. Ausas, Fabricio S. Sousa, Felipe Pereira, Gustavo C. Buscaglia

Abstract: 

A three dimensional parallel implementation of Multiscale Mixed Methods based on non-overlapping domain decomposition techniques is proposed for multi-core computers and its computational performance is assessed by means of numerical experiments. As a prototypical method, from which many others can be derived, the Multiscale Robin Coupled Method is chosen and its implementation explained in detail. Numerical results for problems ranging from millions up to more than 2 billion computational cells in highly heterogeneous anisotropic rock formations based on the SPE10 benchmark are shown. The proposed implementation relies on direct solvers for both local problems and the interface coupling system. We find good weak and strong scalalability as compared against a state-of-the-art global fine grid solver based on Algebraic Multigrid preconditioning in single and two-phase flow problems.

DOI Open Access | ArXiv Preprint | PDF

Applied Mathematics and Computation, v. 421, p. 126908, 2022.

Franciane F. Rocha, Het Mankad, Fabricio S. Sousa, Felipe Pereira  

Abstract: 

In this work we formulate and test a new procedure, the Multiscale Perturbation Method for Two-Phase Flows (MPM-2P), for the fast, accurate and naturally parallelizable numerical solution of two-phase, incompressible, immiscible displacement in porous media approximated by an operator splitting method. The proposed procedure is based on domain decomposition and combines the Multiscale Perturbation Method (MPM) [Ali, et al., Appl. Math. and Comput., 387 (2020) pp. 125023] with the Multiscale Robin Coupled Method (MRCM) [Guiraldello, et al., J. Comput. Phys., 355 (2018) pp. 1-21]. When an update of the velocity field is called for by the operator splitting algorithm, the MPM-2P may provide, depending on the magnitude of a dimensionless algorithmic parameter, an accurate and computationally inexpensive approximation for the velocity field by reusing previously computed multiscale basis functions. Thus, a full update of all multiscale basis functions required by the MRCM for the construction of a new velocity field is avoided. There are two main steps in the formulation of the MPM-2P. Initially, for each subdomain one local boundary value problem with trivial Robin boundary conditions is solved (instead of a full set of multiscale basis functions, that would be required by the MRCM). Then, the solution of an inexpensive interface problem provides the velocity field on the skeleton of the decomposition of the domain. The resulting approximation for the velocity field is obtained by downscaling. We consider challenging two-phase flow problems, with high-contrast permeability fields and water-oil finger growth in homogeneous media. Our numerical experiments show that the use of the MPM-2P gives exceptional speed-up - almost 90% of reduction in computational cost - of two-phase flow simulations. Hundreds of MRCM solutions can be replaced by inexpensive MPM-2P solutions, and water breakthrough can be simulated with very few updates of the MRCM set of multiscale basis functions.

DOI | PDF

Computer Methods in Applied Mechanics and Engineering, v. 385, p. 114035, 2021.

Franciane F. Rocha, Fabricio S. Sousa, Roberto F. Ausas, Felipe Pereira, Gustavo C. Buscaglia  

Abstract: 

It is well known that domain-decomposition-based multiscale mixed methods rely on interface spaces, defined on the skeleton of the decomposition, to connect the solution among the non-overlapping subdomains. Usual spaces, such as polynomial-based ones, cannot properly represent high-contrast channelized features such as fractures (high permeability) and barriers (low permeability) for flows in heterogeneous porous media. We propose here new interface spaces, which are based on physics, to deal with permeability fields in the simultaneous presence of fractures and barriers, accommodated respectively, by the pressure and flux spaces. Existing multiscale methods based on mixed formulations can take advantage of the proposed interface spaces, however, in order to present and test our results, we use the newly developed Multiscale Robin Coupled Method (MRCM) (Guiraldello et al., 2018), which generalizes most well-known multiscale mixed methods, and allows for the independent choice of the pressure and flux interface spaces. An adaptive version of the MRCM (Rocha et al., 2020) is considered that automatically selects the physics-based pressure space for fractured structures and the physics-based flux space for regions with barriers, resulting in a procedure with improved accuracy. The features of the proposed approach are investigated through several numerical simulations of single-phase and two-phase flows, in different heterogeneous porous media. The adaptive MRCM combined with the interface spaces based on physics provides promising results for challenging problems with the simultaneous presence of fractures and barriers.

DOI | PDF

Applied Mathematical Modelling, v. 91, p. 1100-1116, 2021.

Stevens Paz, Alfredo Jaramillo, Rafael T. Guiraldello, Roberto F. Ausas, Fabricio S. Sousa, Felipe Pereira, Gustavo C. Buscaglia  

Abstract:

Two-phase flows in oil reservoirs can be modeled by a coupled system of elliptic and hyperbolic partial differential equations. The transport velocity of the multiphase fluid system is related to the pressure through Darcy’s law and it is coupled to a conservation law for the saturation variable of one of the phases. A time step of the classical IMPES (IMplicit Pressure Explicit Saturation) method consists of first solving the elliptic problem for pressure and Darcy velocity, and then updating the saturation with an explicit numerical scheme for conservation laws. This method is very computationally costly, since the time-consuming elliptic solver must be invoked at time intervals defined by the stability limit of the hyperbolic solver. A popular variant is not to update the velocity at all hyperbolic time steps, but to skip a fixed number C of them, with C determined by the user. In this work we propose a more accurate and systematic procedure for time stepping in IMPES codes. The velocity is updated at all transport time steps, though the elliptic solver is only invoked every C steps. In the time steps at which the elliptic problem is not solved, the velocity is extrapolated from previously computed values with polynomials of high degree. Further, we introduce an error estimator that allows for the number C to be adaptively determined without user intervention. The algorithm was tested in several relevant benchmark problems. This allowed for the optimization of its parameters and comparisons with previous variants. The results show that the proposed algorithm is very stable, reliable and time-cost effective. It is also easily implemented in pre-existent IMPES codes.

DOI | PDF 

Journal of Computational Physics, v. 409, p. 109316, 2020.

Franciane F. Rocha, Fabricio S. Sousa, Roberto F. Ausas, Gustavo C. Buscaglia, Felipe Pereira  

Abstract: 

The Multiscale Robin Coupled Method (MRCM) is a domain decomposition procedure that has been developed to efficiently approximate velocity and pressure fields for single-phase flows in highly heterogeneous porous media. It generalizes other well-established multiscale domain decomposition mixed methods and it adds great flexibility to the choice of interface spaces as well as in the boundary conditions for subdomain coupling. We investigate the approximation of two phase flows in porous media using the MRCM to compute velocity fields. We consider an operator splitting strategy, where a scalar conservation law for the saturation of one of the phases and the velocity field are updated sequentially. We find that the coupling of nearest neighbor subdomains through the imposition of a continuous pressure (respectively, normal fluxes) is the best strategy to approximate two-phase flows in the presence of high (resp., low) permeability channels (resp., regions). A new adaptivity strategy for setting an algorithmic parameter of the MRCM, that controls the relative importance of Dirichlet and Neumann boundary conditions in the coupling of subdomains, is proposed and tested in challenging, high-contrast permeability fields. Our numerical simulations of two-phase flows show that by switching between existing multiscale procedures we can observe unprecedented accuracy, in that we produce better solutions for problems with high-contrast permeability coefficients when compared to solutions obtained with some standard multiscale mixed methods.

DOI | PDF

Computational Geosciences, v. 24, p. 1141–1161, 2020.

Rafael T. Guiraldello, Roberto F. Ausas, Fabricio S. Sousa, Felipe Pereira, Gustavo C. Buscaglia

Abstract: 

We propose two postprocessing procedures (Patch method and Stitch method) to recover local conservation of velocity fields produced by multiscale approximations that are only conservative in coarse scales. These procedures operate on small overlapping regions and are designed to be implemented in parallel, which makes them relatively inexpensive. We investigate the applicability of such methods when tested on single-phase flow problems using the Multiscale Robin Coupled Method (MRCM) in highly heterogeneous permeability fields for modeling the contaminant transport in the subsurface. Numerical simulations are presented aiming to illustrate and compare the performance of the new methods with a standard procedure, the Mean method, in terms of accuracy in contaminant concentration. We show that for a collection of permeability fields taken as log-normal fields, the new postprocessing procedures provide similar or better accuracy than the Mean method. Then, we turn our attention to flows in high-contrast channelized porous formations, where the new methods robustly yield more accurate results and should thus be favored.

DOI | PDF

Applied Mathematics and Computation, v. 387, p. 125023, 2020.

Alsadig Ali, Het Mankad, Felipe Pereira, Fabricio S. Sousa

Abstract: 

In the numerical solution of elliptic equations, multiscale methods typically involve two steps: the solution of families of local solutions or multiscale basis functions (an embarrassingly parallel task) associated with subdomains of a domain decomposition of the original domain, followed by the solution of a global problem. In the solution of multiphase flow problems approximated by an operator splitting method one has to solve an elliptic equation every time step of a simulation, that would require that all multiscale basis functions be recomputed. In this work, we focus on the development of a novel method that replaces a full update of local solutions by reusing multiscale basis functions that are computed at an earlier time of a simulation. The procedure is based on classical perturbation theory. It can take advantage of both an offline stage (where multiscale basis functions are computed at the initial time of a simulation) as well as of a good initial guess for velocity and pressure. The formulation of the method is carefully explained and several numerical studies are presented and discussed. They provide an indication that the proposed procedure can be of value in speeding-up the solution of multiphase flow problems by multiscale methods.

DOI | PDF

Mathematics and Computers in Simulation, v. 164, p. 103-119, 2019.

Rafael T. Guiraldello, Roberto F. Ausas, Fabricio S. Sousa, Felipe Pereira, Gustavo C. Buscaglia

Abstract: 

The Multiscale Robin Coupled Method (MRCM) is a recent multiscale numerical method based on a non-overlapping domain decomposition procedure. One of its hallmarks is that the MRCM allows for the independent definition of interface spaces for pressure and flux over the skeleton of the decomposition. The accuracy of the MRCM depends on the choice of these interface spaces, as well as on an algorithmic parameter β in the Robin interface conditions imposed at the subdomain boundaries. This work presents an extensive numerical assessment of the MRCM in both of these aspects. Two types of interface spaces are implemented: usual piecewise polynomial spaces and informed spaces, the latter obtained from sets of snapshots by dimensionality reduction. Different distributions of the unknowns between pressure and flux are explored. Two nondimensionalizations of β are tested. The assessment is conducted on realistic, high contrast, channelized permeability fields from a SPE benchmark database. The results show that β, suitably nondimensionalized, can be fixed to unity to avoid any indeterminacy in the method. Further, with both types of spaces it is observed that a balanced distribution of the interface unknowns between pressure and flux renders the MRCM quite attractive both in accuracy and in computational cost, competitive with other multiscale methods from the literature.

DOI | PDF

Journal of Computational Physics, v. 396, p. 848-866, 2019.

Fabricio S. Sousa, Camila F. Lages, Jonas L. Ansoni, Antonio Castelo, Adenilso Simao 

Abstract: 

Tree-based mesh grids bring the advantage of using fast cartesian discretizations, such as finite differences, and the flexibility and accuracy of local mesh refinement. The main challenge is how to adapt the discretization stencil near the interfaces between grid elements of different sizes, which is usually solved by local high-order geometrical interpolations. These interpolations depend on the distribution of cells in the vicinity of the point of interest, hence they are site-specific and can become quite complex in three-dimensional simulations, specially when dealing with staggered unknown arrangements. Most methods usually avoid this by limiting the mesh configuration (usually to graded quadtree/octree grids), reducing the number of cases to be treated locally. In this work, we propose a robust method based on a moving least squares meshless interpolation technique, which is employed to compute the weights of the finite difference approximation in a given hierarchical grid, allowing for complex mesh configurations, still keeping the overall order of accuracy of the resulting method. Numerical convergence tests and application to fluid flow simulations are performed to illustrate the flexibility, robustness and accuracy of this new approach.

DOI | PDF

Journal of Computational Physics, v. 355, p. 1-21, 2018.

Rafael T. Guiraldello, Roberto F. Ausas, Fabricio S. Sousa, Felipe Pereira, Gustavo C. Buscaglia

Abstract: 

A multiscale mixed method aiming at the accurate approximation of velocity and pressure fields in heterogeneous porous media is proposed. The procedure is based on a new domain decomposition method in which the local problems are subject to Robin boundary conditions. The domain decomposition procedure is defined in terms of two independent spaces on the skeleton of the decomposition, corresponding to interface pressures and fluxes, that can be chosen with great flexibility to accommodate local features of the underlying permeability fields. The well-posedness of the new domain decomposition procedure is established and its connection with the method of Douglas et al. (1993) [12], is identified, also allowing us to reinterpret the known procedure as an optimized Schwarz (or Two-Lagrange-Multiplier) method. The multiscale property of the new domain decomposition method is indicated, and its relation with the Multiscale Mortar Mixed Finite Element Method (MMMFEM) and the Multiscale Hybrid-Mixed (MHM) Finite Element Method is discussed. Numerical simulations are presented aiming at illustrating several features of the new method. Initially we illustrate the possibility of switching from MMMFEM to MHM by suitably varying the Robin condition parameter in the new multiscale method. Then we turn our attention to realistic flows in high-contrast, channelized porous formations. We show that for a range of values of the Robin condition parameter our method provides better approximations for pressure and velocity than those computed with either the MMMFEM and the MHM. This is an indication that our method has the potential to produce more accurate velocity fields in the presence of rough, realistic permeability fields of petroleum reservoirs.

DOI | PDF

Journal of the Brazilian Society of Mechanical Sciences and Engineering, v. 40, 417, 2018.

Abstract: 

This paper describes in detail a numerical scheme to predict complex turbulent flows using a recent model based on temporal large-eddy simulations (TLES). To solve the equations a second-order finite volume numerical method coupled with a second-order time integration scheme is used. The numerical scheme is validated and then applied to present new results concerning the prediction of the complex turbulent flow in a cubic lid-driven cavity, at Reynolds numbers Re=12,000 and Re=18,000. The results obtained with the TLES are compared with direct numerical simulations and experimental data for the mean velocity flow field and for the Reynolds stresses, showing to be very attractive when compared to large-eddy simulations.

DOI | Link