Multi-frontal Direct Solver Research Articles

Frequency-Domain Reverse-Time Migration with Analytic Green's Function for the Seismic Imaging of Shallow Water Column Structures in the Arctic Ocean.

Seismic oceanography can provide a two- or three-dimensional view of the water column thermocline structure at a vertical and horizontal resolution from the multi-channel seismic dataset. Several seismic imaging methods and techniques for seismic oceanography have been presented in previous research. In this study, we suggest a new formulation of the frequency-domain reverse-time migration method for seismic oceanography based on the analytic Green's function. For imaging thermocline structures in the water column from the seismic data, our proposed seismic reverse-time migration method uses the analytic Green's function for numerically calculating the forward- and backward-modeled wavefield rather than the wave propagation modeling in the conventional algorithm. The frequency-domain reverse-time migration with analytic Green's function does not require significant computational memory, resources, or a multifrontal direct solver to calculate the migration seismic images as like conventional reverse-time migration. The analytic Green's function in our reverse-time method makes it possible to provide a high-resolution seismic water column image with a meter-scale grid size, consisting of full-band frequency components for a modest cost and in a low-memory environment for computation. Our method was applied to multi-channel seismic data acquired in the Arctic Ocean and successfully constructed water column seismic images containing the oceanographic reflections caused by thermocline structures of the water mass. From the numerical test, we note that the oceanographic reflections of the migrated seismic images reflected the distribution of Arctic waters in a shallow depth and showed good correspondence with the anomalies of measured temperatures and calculated reflection coefficients from each XCDT profile. Our proposed method has been verified for field data application and accuracy of imaging performance.

Is 3D frequency-domain FWI of full-azimuth/long-offset OBN data feasible? The Gorgon data FWI case study

Frequency-domain full-waveform inversion (FWI) is potentially amenable to efficient processing of full-azimuth long-offset stationary-recording seabed acquisition carried out with a sparse layout of ocean-bottom nodes (OBNs) and broadband sources because the inversion can be performed with a few discrete frequencies. However, computing the solution of the forward (boundary-value) problem efficiently in the frequency domain with linear algebra solvers remains a challenge for large computational domains involving tens to hundreds of millions of parameters. We illustrate the feasibility of 3D frequency-domain FWI with a subset of the 2015/2016 Gorgon OBN data set in the North West Shelf, Australia. We solve the forward problem with the massively parallel multifrontal direct solver MUMPS, which includes four key features to reach high computational efficiency: an efficient parallelism combining message-passing interface and multithreading, block low-rank compression, mixed-precision arithmetic, and efficient processing of sparse sources. The Gorgon subdata set involves 650 OBNs that are processed as reciprocal sources and 400,000 sources. Monoparameter FWI for vertical wavespeed is performed in the viscoacoustic vertically transverse isotropic approximation with a classical frequency continuation approach proceeding from a starting frequency of 1.7 Hz to a final frequency of 13 Hz. The target covers an area ranging from 260 km2 (frequency ≥ 8.5 Hz) to 705 km2 (frequency ≤ 8.5 Hz) for a maximum depth of 8 km. Compared to the starting model, FWI dramatically improves the reconstruction of the bounding faults of the Gorgon horst at reservoir depths as well as several intrahorst faults and several horizons of the Mungaroo Formation down to a depth of 7 km. Seismic modeling reveals a good kinematic agreement between recorded and simulated data, but amplitude mismatches between the recorded and simulated reflection from the reservoir suggest elastic effects. Therefore, future works involve multiparameter reconstruction for density and attenuation before considering elastic FWI from hydrophone and geophone data.

2.5-Dimensional modeling of EM logging-while-drilling tool in anisotropic medium on a Lebedev grid

A 2.5D finite-difference (FD) algorithm for the modeling of the electromagnetic (EM) logging-while-drilling (LWD) tool in anisotropic media is presented. The FD algorithm is based on the Lebedev grid, which allows for the discretization of the frequency-domain Maxwell's equations in the anisotropic media in 2.5D scenarios without interpolation. This leads to a system of linear equations that is solved using the multifrontal direct solver which enables the simulation of multi-sources at nearly the cost of simulating a single source for each frequency. In addition, near-optimal quadrature derived from an optimized integration path in the complex plane is employed to implement the fast inverse Fourier Transform (IFT). The algorithm is then validated by both analytic and 3D solutions. Numerical results show that two Lebedev subgrid sets are sufficient for TI medium, which is common in geosteering environments. The number of quadrature points is greatly reduced by using the near-optimal quadrature method.

Refined isogeometric analysis of quadratic eigenvalue problems

Certain applications that analyze damping effects require the solution of quadratic eigenvalue problems (QEPs). We use refined isogeometric analysis (rIGA) to solve quadratic eigenproblems. rIGA discretization, while conserving desirable properties of maximum-continuity isogeometric analysis (IGA), reduces the interconnection between degrees of freedom by adding low-continuity basis functions. This connectivity reduction in rIGA’s algebraic system results in faster matrix LU factorizations when using multifrontal direct solvers. We compare computational costs of rIGA versus those of IGA when employing Krylov eigensolvers to solve quadratic eigenproblems arising in 2D vector-valued multifield problems. For large problem sizes, the eigencomputation cost is governed by the cost of LU factorization, followed by costs of several matrix–vector and vector–vector multiplications, which correspond to Krylov projections. We minimize the computational cost by introducing C0 and C1 separators at specific element interfaces for our rIGA generalizations of the curl-conforming Nédélec and divergence-conforming Raviart–Thomas finite elements. Let p be the polynomial degree of basis functions; the LU factorization is up to O((p−1)2) times faster when using rIGA compared to IGA in the asymptotic regime. Thus, rIGA theoretically improves the total eigencomputation cost by O((p−1)2) for sufficiently large problem sizes. Yet, in practical cases of moderate-size eigenproblems, the improvement rate deteriorates as the number of computed eigenvalues increases because of multiple matrix–vector and vector–vector operations. Our numerical tests show that rIGA accelerates the solution of quadratic eigensystems by O(p−1) for moderately sized problems when we seek to compute a reasonable number of eigenvalues.

Accurate 3D frequency-domain seismic wave modeling with the wavelength-adaptive 27-point finite-difference stencil: A tool for full-waveform inversion

Efficient frequency-domain full-waveform inversion (FWI) of long-offset node data can be performed with a few frequencies. The seismic response of these frequencies can be computed with compact finite-difference stencils on regular Cartesian grid with direct or hybrid direct/iterative methods. Compactness, which is necessary to mitigate the fill-in induced by the lower-upper factorization, is implemented with second-order stencils, whereas accuracy is achieved by building a consistent mass matrix and a compound stiffness matrix with optimal weights so that the stencil covers the eight cells surrounding the central point, leading to a 27-point stencil. Classical approaches estimate constant weights by jointly minimizing numerical dispersion in homogeneous media for several numbers of grid points per wavelength ([Formula: see text]). Then, the impedance matrix is built by using these constant weights at each central point covering the heterogeneous subsurface model, leading to nonuniform wavefield accuracy. Instead, we estimate [Formula: see text]-dependent weights once and for all by minimizing dispersion separately for each value of [Formula: see text] in a range found in FWI applications. Then, we build the impedance matrix without computational overhead by selecting at each central point the weights corresponding to the local [Formula: see text], hence leading to a wavelength-adaptive stencil. This separate approach with adaptive weights makes wavefield accuracy uniform. We benchmark wavefields, which are computed in several 3D large-scale subsurface models with a sparse multifrontal direct solver and the nonadaptive/adaptive stencils, against analytical solutions when available and the highly accurate discretization-free convergent Born series method. Each benchmark reveals the higher accuracy of the adaptive stencil relative to the nonadaptive one. In the presence of sharp contrasts, the adaptive stencil also is more accurate than a classical [Formula: see text] finite-different time-domain staggered-grid stencil. The relevance of the adaptive stencil is finally illustrated with a 3D FWI case study in the 3.5–13 Hz frequency band.

DGIRM: Discontinuous Galerkin based isogeometric residual minimization for the Stokes problem

In this paper, we introduce a stable isogeometric analysis discretization of the Stokes system of equations. We use this standard constrained problem to demonstrate the flexibility and robustness of the residual minimization method on dual stable norms [16], which unlocks the extraordinary approximation power of isogeometric analysis [44]. That is, we introduce an isogeometric residual minimization method (IRM) for the Stokes equations, which minimizes the residual in a dual discontinuous Galerkin norm; thus we use the acronym DGiRM. Following Calo et al. [16], we start from an inf-sup stable discontinuous Galerkin (DG) formulation to approximate in a highly continuous trial space that minimizes the dual norm of the residual in a discontinuous test space. We demonstrate the performance and robustness of the methodology considering a manufactured solution and the well-known lid-driven cavity flow problem. First, we use a multi-frontal direct solver, and, using the Pareto front, compare the resulting numerical accuracy and the computational cost expressed by the number of floating-point operations performed by the direct solver algorithm. Second, we use an iterative solver. We measure the number of iterations required when increasing the mesh size and how the configuration of spaces affect the resulting accuracy. This paper is an extension of the paper A Stable Discontinuous Galerkin Based Isogeometric Residual Minimization for the Stokes Problem by M. Łoś, et al. (2020) published in Lecture Notes in Computer Science. In this extended version, we deepen in the mathematical aspects of the DGiRM, include the iterative solver algorithm and implementation, and discuss the influence of different discretization spaces on the iterative solver's convergence.

Computational Costs of Multi-Frontal Direct Solvers with Analysis-Suitable T-Splines

In this paper, we consider the computational cost of a multi-frontal direct solver used for the factorization of matrices resulting from a discretization of isogeometric analysis with T-splines, and analysis-suitable T-splines. We start from model projection or model heat transfer problems discretized over two-dimensional meshes, either uniformly refined or refined towards a point or an edge. These grids preserve several symmetries and they are the building blocks of more complicated grids constructed during adaptive isotropic refinement procedures. A large class of computational problems construct meshes refined towards point or edge singularities. We propose an ordering that permutes the matrix in a way that the computational cost of a multi-frontal solver executed on adaptive grids is linear. We show that analysis-suitable T-splines with our ordering, besides having other well-known advantages, also significantly reduce the computational cost of factorization with the multi-frontal direct solver. Namely, the factorization with N T-splines of order p over meshes refined to a point has a linear O(Np4) cost, and the factorization with T-splines on meshes refined to an edge has O(N2pp2) cost. We compare the execution time of the multi-frontal solver with our ordering to the Approximate Minimum Degree (AMD) and Cuthill–McKee orderings available in Octave.

Real-Time Monitoring and Detection of Emissions in Vehicles

This research aims to optimize the weight of already existing fire extinguishers containing Carbon Dioxide as per BIS standards IS 2190: 2010 and IS 7285-1: 2004. In the paper, an optimized and new design is proposed to minimize the weight of the pressure vessel making it easier to store in smaller compartments and also allowing faster hand-held use during an emergency. The weight has been reduced substantially by using proposed materials and FEA simulations were performed on the model to check for failures under the set test conditions with the aim to minimize the volume of metal used that is able to withstand the hydrostatic load of pressure gas. ANSYS Mechanical 18.2 was used for performing 3d pressure analysis while the model was designed on Solidworks 2016. A Multi Frontal Direct Solver was preferred to reach the solutions in every case. The analysis results of all the materials were compared and finally, a material was selected that was tough enough to clear the safety standards and extremely light so that even old people could use it without any issues. A simple shell of Titanium Grade-5 Alloy was found to be the proper and optimal choice to obtain good mechanical properties and lightweight.

Design and Analysis of Lightweight Fire Extinguishers using Alternately Available Metals

This research aims to optimize the weight of already existing fire extinguishers containing Carbon Dioxide as per BIS standards IS 2190: 2010 and IS 7285-1: 2004. In the paper, an optimized and new design is proposed to minimize the weight of the pressure vessel making it easier to store in smaller compartments and also allowing faster hand-held use during an emergency. The weight has been reduced substantially by using proposed materials and FEA simulations were performed on the model to check for failures under the set test conditions with the aim to minimize the volume of metal used that is able to withstand the hydrostatic load of pressure gas. ANSYS Mechanical 18.2 was used for performing 3d pressure analysis while the model was designed on Solidworks 2016. A Multi Frontal Direct Solver was preferred to reach the solutions in every case. The analysis results of all the materials were compared and finally, a material was selected that was tough enough to clear the safety standards and extremely light so that even old people could use it without any issues. A simple shell of Titanium Grade-5 Alloy was found to be the proper and optimal choice to obtain good mechanical properties and lightweight.

3D frequency-domain elastic wave modeling with the spectral element method using a massively parallel direct solver

ABSTRACTComplex topography, the free-surface boundary condition, and anelastic properties of media should be accounted for in the frame of onshore geophysical prospecting imaging, such as full-waveform inversion (FWI). In this context, an accurate and efficient forward-modeling engine is mandatory. We have performed 3D frequency-domain anisotropic elastic wave modeling by using the highly accurate spectral element method and a sparse multifrontal direct solver. An efficient approach similar to computing the matrix-vector product in the time domain is used to build the matrix. We validate the numerical results by comparing with analytical solutions. A parallel direct solver, the sparse direct multifrontal massively parallel solver (MUMPS), is used to solve the linear system. We find that a hybrid implementation of message passing interface and open multiprocessing is more efficient in flops and memory cost. The influence of the deformed mesh, free-surface boundary condition, and heterogeneity of media on MUMPS performance is negligible. Complexity analysis suggests that the memory complexity of MUMPS agrees with the theoretical order [Formula: see text] (or [Formula: see text] with an efficient matrix reordering method) for an [Formula: see text] grid when nontrivial topography is considered. With the available resources, we conduct a moderate scale modeling with a subset of the SEAM Phase II Foothills model, where 60 wavelengths in the [Formula: see text]-axis are propagated. Computing one gradient of FWI based on this model using the frequency-domain modeling is shown to require similar or fewer computational resources than what would be required for a time-domain solver, depending on the number of sources, while larger memory is necessary. An estimation of the increasing trend indicates that approximately 20 Tb of memory would be required for a [Formula: see text] wavelength modeling. The limit of MUMPS scalability hinders the application to larger scale applications.

HYPERGRAPH GRAMMAR BASED MULTI-THREAD MULTI-FRONTAL DIRECT SOLVER WITH GALOIS SCHEDULER

In this paper we analyze two dimensional grids with point and edge singularities in order to develop an efficient graph grammar based multi-frontal direct solver algorithm. We express these grids by hypergraph models. For these meshes we define a sequence of graph grammar productions expressing the construction of frontal matrices, elimination of fully assembled nodes, merging of resulting Schur complements, and repeating the process of elimination and merging until a single frontal matrix remains. The dependency relation between graph grammar productions is analyzed, and the dependency graph is plot, which is equivalent to the elimination tree of the multi-frontal solver algorithm. We utilize classical multi-frontal solver algorithm, and the graph grammar productions allows us to construct an efficient elimination tree, based on the graph representation of the computational mesh, and not the global matrix itself. The graph grammar productions are assigned to nodes of the dependency graph, and they are implemented as tasks in the GALOIS system and scheduled according to the developed dependency graph over the shared memory parallel machine. We show that our graph grammar based solver outperforms parallel MUMPS solver.

Large-scale 3-D EM modelling with a Block Low-Rank multifrontal direct solver

We put forward the idea of using a Block Low-Rank (BLR) multifrontal direct solver to efficiently solve the linear systems of equations arising from a finite-difference discretization of the frequency-domain Maxwell equations for 3-D electromagnetic (EM) problems. The solver uses a low-rank representation for the off-diagonal blocks of the intermediate dense matrices arising in the multifrontal method to reduce the computational load. A numerical threshold, the so-called BLR threshold, controlling the accuracy of low-rank representations was optimized by balancing errors in the computed EM fields against savings in floating point operations (flops). Simulations were carried out over large-scale 3-D resistivity models representing typical scenarios for marine controlled-source EM surveys, and in particular the SEG SEAM model which contains an irregular salt body. The flop count, size of factor matrices and elapsed run time for matrix factorization are reduced dramatically by using BLR representations and can go down to, respectively, 10, 30 and 40 per cent of their full-rank values for our largest system with N = 20.6 million unknowns. The reductions are almost independent of the number of MPI tasks and threads at least up to 90 × 10 = 900 cores. The BLR savings increase for larger systems, which reduces the factorization flop complexity from O(N2) for the full-rank solver to O(Nm) with m = 1.4–1.6. The BLR savings are significantly larger for deep-water environments that exclude the highly resistive air layer from the computational domain. A study in a scenario where simulations are required at multiple source locations shows that the BLR solver can become competitive in comparison to iterative solvers as an engine for 3-D controlled-source electromagnetic Gauss–Newton inversion that requires forward modelling for a few thousand right-hand sides.

APPLICATIONS OF ALTERNATING DIRECTION SOLVER FOR SIMULATIONS OF TIME-DEPENDENT PROBLEMS

This paper deals with application of Alternating Direction solver (ADS) to nonstationarylinear elasticity problem solved with isogeometric FEM. Employingtensor product B-spline basis in isogeometric analysis under some restrictionsleads to system of linear equations with matrix possessing tensor product structure.Alternating Direction Implicit algorithm is a direct method that exploitsthis structure to solve the system in O (N ), where N is a number of degreesof freedom (basis functions). This is asymptotically faster than state-of-theartgeneral purpose multi-frontal direct solvers. In this paper we also presentthe complexity analysis of ADS incorporating dependence on order of B-splinebasis.

Fast 3D frequency-domain full-waveform inversion with a parallel block low-rank multifrontal direct solver: Application to OBC data from the North Sea

Wide-azimuth long-offset ocean bottom cable (OBC)/ocean bottom node surveys provide a suitable framework to perform computationally efficient frequency-domain full-waveform inversion (FWI) with a few discrete frequencies. Frequency-domain seismic modeling is performed efficiently with moderate computational resources for a large number of sources with a sparse multifrontal direct solver (Gauss-elimination techniques for sparse matrices). Approximate solutions of the time-harmonic wave equation are computed using a block low-rank (BLR) approximation, leading to a significant reduction in the operation count and in the volume of communication during the lower upper (LU) factorization as well as offering great potential for reduction in the memory demand. Moreover, the sparsity of the seismic source vectors is exploited to speed up the forward elimination step during the computation of the monochromatic wavefields. The relevance and the computational efficiency of the frequency-domain FWI performed in the viscoacoustic vertical transverse isotropic (VTI) approximation was tested with a real 3D OBC case study from the North Sea. The FWI subsurface models indicate a dramatic resolution improvement relative to the initial model built by reflection traveltime tomography. The amplitude errors introduced in the modeled wavefields by the BLR approximation for different low-rank thresholds have a negligible footprint in the FWI results. With respect to a standard multifrontal sparse direct factorization, and without compromise of the accuracy of the imaging, the BLR approximation can bring a reduction of the LU factor size by a factor of up to three. This reduction is not yet exploited to reduce the effective memory usage (ongoing work). The flop reduction can be larger than a factor of 10 and can bring a factor of time reduction of around three. Moreover, this reduction factor tends to increase with frequency, namely with the matrix size. Frequency-domain viscoacoustic VTI FWI can be viewed as an efficient tool to build an initial model for elastic FWI of 4C OBC data.

Hybrid Direct and Iterative Solver with Library of Multi-criteria Optimal Orderings for h Adaptive Finite Element Method Computations

In this paper we present a multi-criteria optimization of element partition trees and resulting orderings for multi-frontal solver algorithms executed for two dimensional h adaptive finite element method. In particular, the problem of optimal ordering of elimination of rows in the sparse matrices resulting from adaptive finite element method computations is reduced to the problem of finding of optimal element partition trees. Given a two dimensional h refined mesh, we find all optimal element partition trees by using the dynamic programming approach. An element partition tree defines a prescribed order of elimination of degrees of freedom over the mesh. We utilize three different metrics to estimate the quality of the element partition tree. As the first criterion we consider the number of floating point operations(FLOPs) performed by the multi-frontal solver.As the second criterion we consider the number of memory transfers (MEMOPS) performed by the multi-frontal solver algorithm. As the third criterion we consider memory usage (NONZEROS) of the multi-frontal direct solver. We show the optimization results for FLOPs vs MEMOPS as well as for the execution time estimated as FLOPs+100*MEMOPS vs NONZEROS. We obtain Pareto fronts with multiple optimal trees, for each mesh, and for each refinement level. We generate a library of optimal elimination trees for small grids with local singularities. We also propose an algorithm that for a given large mesh with identified local sub-grids, each one with local singularity. We compute Schur complements over the sub-grids using the optimal trees from the library, and we submit the sequence of Schur complements into the iterative solver ILUPCG.

Bisection weighted by element size ordering algorithm for multi-frontal solver executed over 3D h refined grids

In this paper we present an algorithm for generation of ordering over 3D grids h refined towards singularities. The ordering controls the execution of multi-frontal direct solver algorithm on systems of linear equations generated by 3D h adaptive finite element method. Our algorithm uses additional knowledge about the structure of the computational mesh, not available to alternative ordering algorithms. The proposed ordering algorithm outperforms other state-of-the-art orderings available through MUMPS interface, namely nested-dissections, AMD and PORD.

Scalability of Direct Solver for Non-stationary Cahn-Hilliard Simulations with Linearized time Integration Scheme

We study the features of a new mixed integration scheme dedicated to solving the non-stationary variational problems. The scheme is composed of the FEM approximation with respect to the space variable coupled with a 3-leveled time integration scheme with a linearized right-hand side operator. It was applied in solving the Cahn-Hilliard parabolic equation with a nonlinear, fourth-order elliptic part. The second order of the approximation along the time variable was proven. Moreover, the good scalability of the software based on this scheme was confirmed during simulations. We verify the proposed time integration scheme by monitoring the Ginzburg-Landau free energy. The numerical simulations are performed by using a parallel multi-frontal direct solver executed over STAMPEDE Linux cluster. Its scalability was compared to the results of the three direct solvers, including MUMPS, SuperLU and PaSTiX.

Computational cost of isogeometric multi-frontal solvers on parallel distributed memory machines

This paper derives theoretical estimates of the computational cost for isogeometric multi-frontal direct solver executed on parallel distributed memory machines. We show theoretically that for the Cp−1 global continuity of the isogeometric solution, both the computational cost and the communication cost of a direct solver are of order O(log(N)p2) for the one dimensional (1D) case, O(Np2) for the two dimensional (2D) case, and O(N4/3p2) for the three dimensional (3D) case, where N is the number of degrees of freedom and p is the polynomial order of the B-spline basis functions. The theoretical estimates are verified by numerical experiments performed with three parallel multi-frontal direct solvers: MUMPS, PaStiX and SuperLU, available through PETIGA toolkit built on top of PETSc. Numerical results confirm these theoretical estimates both in terms of p and N. For a given problem size, the strong efficiency rapidly decreases as the number of processors increases, becoming about 20% for 256 processors for a 3D example with 1283 unknowns and linear B-splines with C0 global continuity, and 15% for a 3D example with 643 unknowns and quartic B-splines with C3 global continuity. At the same time, one cannot arbitrarily increase the problem size, since the memory required by higher order continuity spaces is large, quickly consuming all the available memory resources even in the parallel distributed memory version. Numerical results also suggest that the use of distributed parallel machines is highly beneficial when solving higher order continuity spaces, although the number of processors that one can efficiently employ is somehow limited.

Volume and neighbors algorithm for finding elimination trees for three dimensional [formula omitted]-adaptive grids

This paper presents an algorithm called “volume & neighbors” for finding elimination trees for multifrontal solver algorithm applied to three dimensional h-adaptive finite element method computations. The algorithm is described in a pseudo-code and explained on exemplary h-refined grids. The algorithm is implemented in three dimensional h-adaptive finite element method code and tested on a sequence of representative grids, namely uniform grid, grid with point singularity, grid with edge singularity and grid with face singularity. The number of floating point operations for the multifrontal solver algorithm working with the elimination trees generated by the volume & neighbors algorithm is compared with the number of floating-point operations resulting from execution of the state-of-the-art multifrontal direct solver MUMPS with state-of-the-art algorithms for constructing elimination trees like nested-dissection from METIS library and approximate minimum degree AMD algorithm as well as PORD algorithm. In all the cases the volume & neighbors algorithm outperforms other state-of-the art algorithms. The only exception is the case of the uniform grid, where the algorithm results in a similar number of FLOPs than nested dissection algorithm.

Dynamic Programming Algorithm for Generation of Optimal Elimination Trees for Multi-frontal Direct Solver Over H-refined Grids

In this paper we present a dynamic programming algorithm for finding optimal elimination trees for computational grids refined towards point or edge singularities. The elimination tree is utilized to guide the multi-frontal direct solver algorithm. Thus, the criterion for the optimization of the elimination tree is the computational cost associated with the multi-frontal solver algorithm executed over such tree. We illustrate the paper with several examples of optimal trees found for grids with point, isotropic edge and anisotropic edge mixed with point singularity. We show the comparison of the execution time of the multi-frontal solver algorithm with results of MUMPS solver with METIS library, implementing the nested dissection algorithm.