key: cord-0560006-uap1vvbm
authors: Grave, Mal'u; Coutinho, Alvaro L. G. A.
title: Adaptive Mesh Refinement and Coarsening for Diffusion-Reaction Epidemiological Models
date: 2020-10-22
journal: nan
DOI: nan
sha: be85815a3b040c5caf93dd7e8861ac99d4e38f80
doc_id: 560006
cord_uid: uap1vvbm

The outbreak of COVID-19 in 2020 has led to a surge in the interest in the mathematical modeling of infectious diseases. Disease transmission may be modeled as compartmental models, in which the population under study is divided into compartments and has assumptions about the nature and time rate of transfer from one compartment to another. Usually, they are composed of a system of ordinary differential equations (ODE's) in time. A class of such models considers the Susceptible, Exposed, Infected, Recovered, and Deceased populations, the SEIRD model. However, these models do not always account for the movement of individuals from one region to another. In this work, we extend the formulation of SEIRD compartmental models to diffusion-reaction systems of partial differential equations to capture the continuous spatio-temporal dynamics of COVID-19. Since the virus spread is not only through diffusion, we introduce a source term to the equation system, representing exposed people who return from travel. We also add the possibility of anisotropic non-homogeneous diffusion. We implement the whole model in texttt{libMesh}, an open finite element library that provides a framework for multiphysics, considering adaptive mesh refinement and coarsening. Therefore, the model can represent several spatial scales, adapting the resolution to the disease dynamics. We verify our model with standard SEIRD models and show several examples highlighting the present model's new capabilities.

Disease transmission may be modeled as compartmental models, in which the population under study is divided into compartments and has assumptions about the nature and time rate of transfer from one compartment to another [14] . These models have been used extensively in biological, ecological, and chemical applications [15, 16, 17] . They allow for an understanding of the processes at work and for predicting the dynamics of the epidemic.

The large majority of the compartmental models are composed of a system of ordinary differential equations (ODE's) in time. Though compartmentalized models are simple to formulate, analyze, and solve numerically, these models do not always account for the movement of individuals from one region to another. Different approaches have been used to introduce spatial variation into such ODE models [18, 19, 20, 11] . The strategies consist of defining regional compartments corresponding to different geographic units and adding coupling terms to the model equations to account for species' movement from unit to unit.

In this work, we use a partial differential equation (PDE) model to capture the continuous spatio-temporal dynamics of COVID-19. PDE models incorporate spatial information more naturally and allow for capture the dynamics across several scales of interest. They have a significant advantage over ODE models, whose ability to describe spatial information is limited by the number of geographic compartments. Indeed, recent research indicates that COVID-19 spreading presents multi-scale features that go from the virus and individual immune system scale to the collective behavior of a whole population [21] . We study a compartmental SEIRD model (susceptible, exposed, infected, recovered, deceased) that incorporates spatial spread through diffusion terms [16, 22, 8, 9, 23] . Adaptive mesh refinement and coarsening [24] can resolve population dynamics from local (street, city) to regional (district, state), providing an accurate spatio-temporal description of the infection spreading. Moreover, diffusion may be properly tuned to account for local natural or social inhomogeneities (e.g., mountains, lakes, highways) describing populations' movements.

However, the main limitation of the diffusion-reaction PDE approach is the definition of the diffusion operator and transmission coefficients, which depend on the population's behavior. Another issue is that the virus spread is not only through diffusion, since people, who may be infected, travel long distances in a short period. Some models relate the mobile geolocation data with the spread of the disease [25, 26] . These issues make the model a highly complex system, which may completely change as the population's behavior changes. Therefore, this work contributes to improving the knowledge of compartmental diffusion-reaction PDE models.

All implementations are done using the libMesh library. As other freely available open-source libraries (deal.II [27] , FEniCS [28] , GRINS [29] , MOOSE [30] , etc), libmesh provides a finite element framework that can be used for numerical simulation of partial differential equations in various applied fields in science and engineering. It has already been used in more than 1,000 publications with applications in many different areas. See, for example, recent applications in sediment transport [31] and bubble dynamics [32] . This library is an excellent tool for programming the finite element method and can be used for one-, two-, and three-dimensional steady and transient simulations on serial and parallel platforms. The libmesh library provides native support for adaptive mesh refinement and coarsening, thus providing a natural environment for the present study. The main advantage of this library is the possibility of focusing on implementing the specifics features of the modeling without worrying about adaptivity and code parallelization. Consequently, the effort to build a high performance computing code tends to be minimized. about adaptivity and code parallelization.

The remainder of this work is organized as follows: In section 2, we present the governing equations that describe the dynamics of a virus infection. First, we present a generic spatio-temporal SEIRD model, based on the EPIDEMIC software [33] , used to verify our implementation. We then present a model that better represents the dynamics of COVID-19 infection spread, based on [9, 8] . In section 3, we introduce the Galerkin finite element formulation, the time discretization, and the libMesh implementation. Then, we present the numerical verification of the generic spatio-temporal SEIRD model implementation. We verify our algorithm's capacity to represent a compartmental model [33] and show how the diffusion influences the dynamics. Section 5 presents the numerical results of the spatio-temporal model of COVID-19 infection spread. We perform simulations similar to the ones presented in [8] and show tests to highlight the new modeling capabilities introduced in this work. Finally, the paper ends with a summary of our main findings and the perspectives for the next steps of this research.

The presentation of the governing equations follows the continuum mechanics framework in [8] instead of the more traditional approach found in mathematical and biological references. Consider a system which may be decomposed into N distinct populations: u 1 (x, t), u 2 (x, t), ..., u N (x, t). Let Ω ∈ R 2 be a simply connected domain of interest with boundary ∂Ω = Γ D ∩ Γ N , and [0, T ] a generic time interval. The vector compact representation of the governing equations as a transient nonlinear diffusion-reaction system of equations reads,

We denote the densities of the susceptible, exposed, infected, recovered and deceased populations as s(x, t), e(x, t), i(x, t), r(x, t), and d(x, t). Also, let c(x, t) denote the cumulative number of infected and n(x, t) the sum of the living population; i.e., n(x, t) = s(x, t) + e(x, t) + i(x, t) + r(x, t). We consider u = [s, e, i, r, d] T . The matrices A, B and ν, and the vector f depend on a particular form of the system dynamics. Furthermore, in general, ν = ν(x), that is, diffusion is heterogeneous and anisotropic. Besides the boundary conditions (2), (3), we specify the initial condition

for i = 1, 2, · · · , N .

We first consider a SEIRD model [14] given by the following system of coupled PDEs over Ω × [0, T ]:

where β is transmission rate (days −1 ), α the latent rate (days −1 ), γ the recovery rate (days −1 ), δ the death rate (days −1 ), and ν s , ν e , ν i , ν r are diffusion parameters respectively corresponding to the different population groups (km 2 persons −1 days −1 ). We append to the system of equations homogeneous Neumann boundary conditions, that is, (ν · ∇u) · n = 0.

We can reframe this model in the general form given by equation (1) . Thus, the matrices A, B, ν and the vector f reads, 

This model is based on the EPIDEMIC software 1 , and it is employed to verify our implementation. The system of equations represents that the susceptible population decreases as the exposed population increases. This variation depends on the transmission rate between infected and susceptible. The number of exposed increases because of the transmission rate and decreases when the exposed individuals become infected (after the incubation period). The number of infected increases after the incubation period and decreases depending on the recovery and death rate. The number of deaths depends only on the death rate as the number of recovered depends only on the recovery rate. Finally, the cumulative number of infected depends only on the exposed and the incubation period. The diffusion parameters are included in the model to spread the disease spatially.

Summarizing, this model assumes:

• Movement is proportional to population size; i.e., more movement occurs within heavily populated regions;

• No movement occurs among the deceased population;

• There is a latency period between exposure and the development of symptoms;

• The probability of contagion is inversely proportional to the population size; • The exposed persons will ever develop symptoms;

• Only infected persons are capable of spreading the disease;

• The non-virus mortality rate is not considered in this model;

• New births are not considered in this model.

Note that the EPIDEMIC model's dynamics does not represent the actual COVID19 dynamics since, in the case of COVID19, the exposed population may be asymptomatic and recover without becoming infected and still spread the virus. Thus, a better model would be the one based on [9, 8] .

We begin by making several model assumptions to represent the COVID-19 infection spread adequately [8] :

• Only mortality due the COVID-19 is considered;

• New births are not considered in this model.

• Some portion of exposed persons never develop symptoms, and move directly from the exposed compartment to the recovered compartment (asymptomatic cases);

• Both asymptomatic (exposed) and symptomatic (infected) patients are capable of spreading the disease;

• There is a latency period between exposure and the development of symptoms;

• It is possible that new cases of exposed people appear randomly in the system (exposed people who return from a travel);

• The probability of contagion increases with population size (Allee effect [9] );

• Movement is proportional to population size; i.e., more movement occurs within heavily populated regions;

• No movement occurs among the deceased population; 1 https://americocunhajr.github.io/EPIDEMIC/ [33] Then, the system of equations becomes:

where A characterizes the Allee effect (persons), that takes into account the tendency of outbreaks to cluster around large populations, β i is the transmission rate between symptomatic and susceptible (persons −1 days −1 ), β e is the transmission rate between asymptomatic and susceptible (persons −1 days −1 ), f is a source function that depends on space and time (persons), α is the latent rate (days −1 ), γ e is the recovery rate of the asymptomatic (days −1 ), γ i is the recovery rate of the symptomatic (days −1 ), δ is the death rate (days −1 ), and ν s , ν e , ν i , ν r are the diffusion parameters respectively corresponding to the different population groups (km 2 persons −1 days −1 ). Now, we call exposed who has contact with the virus but remains asymptomatic. However, since the virus is highly transmissible, the exposed population also may transmit the virus. The exposed may recover without any symptoms or may become infected. The infected follow the same logic of the previous SEIRD system (they may recover or die). The main difference in the new SEIRD system is in the exposed population and whom it interacts. The source function f may be defined to represent exposed people who return from travel. Note that β has units (days −1 ) while β i and β e have units (person −1 days −1 ). While equations (5) and (6) divide β by the living population, equations (15), (16) and (17) keep β i and β e constant independent of that.

Therefore, to express this model in the general form given by equation (1), the matrices A, B, ν and the vector f reads,

If we assume that the region of interest is isolated, we prescribe the following homogeneous Neumann boundary conditions,

or simply (ν · ∇u) · n = 0.

The basic reproduction number, R 0 , is defined as the average number of additional infections produced by an infected individual in a wholly susceptible population over the full course of the disease outbreak. In an epidemic situation, the threshold R 0 = 1 is the dividing line between the infection dying out and the onset of an epidemic. R 0 > 1 implies growth of the epidemic, whereas R 0 < 1 implies decay in infectious spread [14] .

The concept of R 0 is well-defined for ODE models. However, its extension to a PDE model is unclear, owing to the influence of diffusion. Viguerie et al. [8] found that a R 0 derived for the ODE version of the PDE model is not consistently reliable to represent the epidemic's dynamic growth adequately. If we do not consider the diffusion, R 0 may be calculated as:

For further details about the R 0 calculation, refer to [34, 8] .

In this section we briefly introduce the Galerkin finite element formulation, the time discretization, and the the libMesh implementation, supporting adaptive mesh refinement and coarsening. Appendices A and B give respectively the resulting finite element matrices for the generic spatio-temporal SEIRD and COVID-19 models.

We introduce a Galerkin finite element variational formulation for space discretization. Without loss of generality, we consider the case of homogeneous Dirichlet and Neumann boundary conditions. Let V u h be a finite dimensional space such that,

in which u h (·, t) is the discrete counterpart of u and w h the weight function. The weak formulation is then: find

Here we define the operation (·, ·) as the standard scalar product in L 2 (Ω).

The SEIRD and COVID-19 models yield stiff systems of equations, making explicit time-marching methods unfeasible. The Backward Euler method is widely applied because of its unconditional numerical stability characteristics. However, it has the disadvantage of being only first-order accurate, which introduces a significant amount of numerical diffusion. Then, we use the second-order Backward Differentiation Formula (BDF2), which, compared to the prevailing Backward Euler method, has significantly better accuracy while retaining unconditional linear stability. The model becomes,

The subscript n + 1 is associated to t = t n+1 and n, and n − 1 to the previous time-steps.

We implement the compartmental epidemiological models in libMesh, a C++ FEM open-source software library for parallel adaptive finite element applications [35] . libMesh also interfaces with external solver packages like PETSc [36] and Trilinos [37] . Recently, libMesh was also coupled with in-situ visualization and data-analysis tools [38, 39] . It provides a finite element framework that can be used for the numerical simulation of partial differential equations on serial and parallel platforms. This library is an excellent tool for programming the finite element method and can be used for one-, two-, and three-dimensional steady and transient simulations. The libMesh library also has native support for adaptive mesh refinement and coarsening (AMR/C).

Multiple scales can be resolved by AMR/C. libMesh supports AMR/C by h-refinement (element subdivision), prefinement (increasing the polynomial approximation order), and hp-refinement, that is, a combination of both [24] . In libMesh, coarsening is supported in the h, p, and hp AMR/C options. In the present work, we restrict ourselves to h-refinement with hanging nodes. The AMR/C procedure uses a local error estimator to drive the refinement and coarsening procedure, considering the error of an element relative to its neighbor elements in the mesh. This error may come from any variable of the system. As it is standard in libMesh, Kelly s error indicator is employed, which uses the H 1 seminorm to estimate the error [40] . Apart from the element interior residual, the flux jumps across the inter-element edges influence the element error. The flux jump of each edge is computed and added to the element error contribution. For both the element residual and flux jump, the values of the desired variables at each node are necessary. Therefore, the error e 2 can be stated as,

where e 2 i is the error of each variable. In this study, we use all population types as variables for the error estimator. After computing the error values, the elements are "flagged" for refining and coarsening regarding their relative error. This is done by a statistical element flagging strategy. It is assumed that the element error e is distributed approximately in a normal probability function. Here, the statistical mean µ s and standard deviation σ s of all errors are calculated. Whether an element is flagged is depending on a refining (r f ) and a coarsening (c f ) fraction. For all errors e < µ s − σ s c f the elements are flagged for coarsening and for all e > µ s + σ s r f the elements are marked for refinement (see Figure 1 ). The refinement level is limited by a maximum h-level (h max ), (see Figure 2) , and the coarsening is done by h-restitution of sub-elements [24] , [41] .

To verify the implementation of the generic spatio-temporal SEIRD model, we have done several tests. For this, we consider a square domain of 1km × 1km centered at (0, 0) for all tests in this section.

In the first test, we do not consider diffusion. We consider a population of 1000 people/km 2 with 1 person/km 2 initially infected in all area of the domain. Then, the initial conditions are: s 0 = 999, e 0 = 0, i 0 = 1, r 0 = 0 and d 0 = 0. This test aims to reproduce a compartmental simulation of the EPIDEMIC software by using the same initial parameters. The results have to be the same in each point of the domain and the same as the EPIDEMIC software. We set α = 0.14286 days −1 , β = 0.25 days −1 , δ = 0.06666 days −1 , γ = 0.1 days −1 and ∆t = 1 day. The mesh has 50 × 50 bilinear quadrilateral elements. Figure 3 shows the comparison of the results, where we can see a very good agreement between both solutions. Figure 4 shows the results over a centralized horizontal line crossing the domain at t=365 days. It is possible to see that the results are the same in all the domain, as expected.

Now, we consider the same parameters of the previous example, but different initial conditions. We consider a population of 1000 people/km 2 in all area of the domain with 1 person/km 2 initially infected only in a circle centered at (0, 0) and radius R = 0.5 km, We assume that ν s = ν e = ν i = ν r = 10 −8 km 2 persons −1 days −1 . Then, the initial conditions are: s 0 = 999, e 0 = 0, i 0 = 1 for R <= 0.5 and i 0 = 0 for R > 0 with R = x 2 + y 2 , r 0 = 0 and d 0 = 0 (see Figure 5 ). We consider adaptive mesh refinement in this example. The original mesh has 50 × 50 bilinear quadrilateral elements, and after the refinement, the smallest element has size 0.005 km. We initially refine the domain in two levels. For the AMR/C procedure, we set h max = 2,r f = 0.95, c f = 0.05. We apply the adaptive mesh refinement every 5 time-steps. Figure 6 shows the results over a centralized horizontal line crossing the domain at t=365 days. Figure 7 shows the infected people at different time-steps. Note that the infected remains actives in other parts of the domain because of the diffusion. It is possible to see the wave effect of the disease spreading. Note that the AMR/C procedure improves spatial resolution on the regions where the infected people are higher. 

In this test, we change the initial population. Instead of a constant value in all domain, we set 1000 people/km 2 at the left/top quadrant, 500 people/km 2 at the right/top quadrant, 250 people/km 2 at the left/bottom quadrant and 750 people/km 2 at the right/bottom quadrant ( Figure 8) . Then, the initial conditions are: s 0 = 999 for x ≤ 0 and y > 0, s 0 = 499 for x > 0 and y > 0, s 0 = 249 for x ≤ 0 and y <=, s 0 = 749 for x > 0 and y > 0, e 0 = 0, i 0 = 1 for R ≤ 0.5 and i 0 = 0 for R > 0 with R = x 2 + y 2 , r 0 = 0 and d 0 = 0. The initial population infected is 1 person/km 2 at the same circled region of the previous test. All other parameters are the same of the previous simulation. Figure 9 shows the infected people ate different time-steps. It is possible to see that the regions with denser populations (more people/km 2 ) are more affected by the disease. Figure 10 shows the total number of deaths after 365 days, and the regions with more people/km 2 have more deaths than the less dense regions. Note also that the AMR/C procedure generates meshes following the model dynamics. 

In this section, we perform some simulations to validate the spatio-temporal model of COVID-19 infection spread.

In this test, we do not consider diffusion. We consider a square domain of 1km × 1km centered at (0, 0) with a population of 1000 people/km 2 , with 1 person/km 2 initially infected and 5 people/km 2 exposed in all area of the domain. Then, the initial conditions are: s 0 = 994, e 0 = 5, i 0 = 1, r 0 = 0 and d 0 = 0. The aim of this test is to reproduce a compartmental simulation presented in [8] by using the same initial parameters. The results has to be the same in each point of the domain and also the same of the ones given in [8] . We set α = 0.125 days −1 , β i = β e = 0.005 days −1 persons −1 , δ = 0.0625 days −1 , γ i = 0.041666667 days −1 and γ e = 0.1666667 days −1 . The mesh has 50 × 50 bilinear quadrilateral elements. Figure 11 shows the comparison of the results, where we can see an excellent agreement. To reproduce a 1D simulation with quadrilateral elements, we fix the element width to 0.0005 and vary its length to find the proper refinement for this case. Therefore, we run a mesh convergence study as well as a time-step convergence study.

For the initial conditions, we set s = s 0 and e = e 0 as follows, (36) Figure 12 shows the initial conditions. We further set i 0 = 0, r 0 = 0, and d 0 = 0. Qualitatively, the initial conditions represent a large population centered around x = 0.35 with no exposed persons and a small population centered around x = 0.75 with some exposed individuals. We also enforce homogeneous Neumann boundary conditions at x = 0 and a zero-population Dirichlet boundary condition at x = 1 for all model compartments. The latter represents a non-populated area at x = 1. 

We compare numerical solutions computed on successively refined uniform grids with mesh size ∆x = 1/50, 1/100, 1/250, 1/500, and 1/1000. The time step is ∆t = 0.25 days. Figure 15 shows the difference in the total population of each compartment of individuals for the different meshes.

A good resolution is found for ∆x = 1/500. It is easy to see this convergence in Figure 16 , where the number of individuals of each compartment is plotted at t = 90 days.

We examine the impact of time-step size ∆t on the numerical approximation of the model solution. We consider the time step sizes ∆t = 1, ∆t = 0.5, ∆t = 0.25, ∆t = 0.125 and ∆t = 0.0625 days. As the results in Section 5.2.1 suggested ∆x = 1/500 is a sufficiently fine spatial discretization, we utilize this mesh resolution here. Figure 17 shows the difference of the total population of each compartment of individuals for the different time-steps.

A good accuracy is found for ∆t = 0.25 days. It is easy to see how the accuracy improves in Figure 18 , where the number of individuals of each compartment is plotted at t = 90 days.

This test is the application of the previous configuration rotated in a two dimensional square with corners at (-1,-1), (1,-1), (1,1) and (-1,1) . The initial population is: 

with R = x 2 + y 2 .

The original mesh has 50 × 50 bilinear quadrilaterals elements and it is refined in two levels at the beginning of the simulation. For the AMR/C procedure, we set h max = 2, r f = 0.95, c f = 0.05. We apply the adaptive mesh refinement every 4 time-steps. The behavior of the transmission has to be similar to the 1D model results, but in a radial configuration. Figures 19 shows the populations at different time steps. Figure 6 shows the results over a centralized horizontal line (or vertical because the axisymmetry) crossing the domain at t=200 days. If we compare Figure 6 with Figure 14 , it is possible to see that the populations follow a similar behavior. In Figure 21 we plot the time history of the total number of individuals. There is a small gain in the total number of individuals (less than 0.1%). 

This test considers anisotropic diffusion in the previous configuration (only in the x direction). Therefore, the populations move spatially only in the x direction. Figure 22 shows the populations at different time-steps. Figure 23 shows the results over a centralized horizontal line crossing the domain, and Figure 24 over a centralized vertical line. By comparing these two figures, it is clear how the diffusion direction influences the behavior of the virus spread. Since there is no movement of infected or exposed people in the y direction, part of the population does not have contact with the virus because there is no chance of the virus to reach them.

In Figure 25 we plot the time history of the total number of individuals. We can see a gain in the total number of individuals of less than 0.1%. The original mesh has 50 × 50 bilinear quadrilaterals elements and it is refined in two levels at the beginning of the simulation. For the AMR/C procedure, we set h max = 2, r f = 0.95, c f = 0.05. We apply the adaptive mesh refinement every 4 time-steps.

The initial population is: Figure 26 shows the initial susceptible population. Note there are not infected or exposed people at the initial time.

We implement a random source of the exposed population that depends on the number of susceptible people. In all time-steps random nodes of the domain receive a certain number of exposed people. It tries to simulate people who travel and suddenly appear in a region carrying the virus. The random source does not add individuals to the population, but change individuals between susceptible and exposed compartments. Of course, this model is simple. Nevertheless, it demonstrates how to handle a random source term in the equations. Figure 27 shows a example of the random exposed number of people that appears in one time-step. In Figure 31 , we plot the time history of the total number of individuals. There is a negligible increase in the total number of individuals (less than 0.1%). 

We developed an extended continuum SEIRD model to represent the dynamics of the COVID-19 virus spread based on the framework proposed in [9] . We validate our code by comparing our results with other simulations. We introduce new test cases to highlight new modeling capabilities. Among the new features added to the base model, there is the addition of a source term, which represents exposed people who return from travel, by changing individuals from the susceptible compartment to the exposed compartment. We also add the possibility of anisotropic non-homogeneous diffusion. Our code is implemented through the libMesh library and supports adaptive mesh refinement and coarsening. Therefore, it can represent several spatial scales, adapting the resolution to the disease dynamics.

Data is essential to define the epidemic spreading parameters, as diffusion and infection rate. We have to study how it would be the best way to represent people who return from travels, addressing questions like defining the probability of a random source appearing in the system, in which area, the population density, among others. Diffusion-reaction models, as the present one, are richer than standard compartmental models. However, they are slower, which hampers their widespread utilization in what-if scenarios, parametric studies, and time-critical situations. Therefore, the development of low-dimensional computational models will leverage the ability of continuous models to perform in real-time scenarios. Projection-based or data-driven model order reduction [42, 43] aims to lower the computational complexity of a given computational model by reducing its dimensionality (or order), can provide this leverage. They can work in conjunction with emerging machine learning methods such as physics informed neural networks [44] . We can foresee a tremendous impact in the mathematical epidemiology field of all these new methods and techniques, enlarging the predictive capabilities and computational efficiency of diffusion-reaction epidemiological models.

In libMesh, we calculate directly the new solution (u n+1 ) instead of the variation (δu). Then, on the left-hand side, we gather the terms containing an unknown, whereas all the other terms are taken to the right-hand-side. The superscript k is from the previous Newton iteration. The terms in black are from the mass matrix, in blue are the nonlinear terms, in red the diffusive terms, and in green the remaining terms from the stiffness matrix. The finite element shape functions are represented by N a , a = 1, · · · , n nnos , where n nnos is the number of nodes in the finite element mesh.

Susceptible (Equation 5): We present the matrix contributions of the system of equations that represents the COVID19 dynamics [9, 8] . We use the BDF2 time discretization method, Newton's method for the nonlinear terms, and we simplify the number of the living population by considering the previous linear solution. For all test cases the nonlinear tolerance for Newton's method is set to 10 −8 and the linear solver tolerance is set to 10 −10 . The linear solver is GMRES with ILU(0) preconditioner.

In libMesh, we calculate directly the new solution (u n+1 ) instead of the variation (δu). Then, on the left-hand side, we gather the terms containing an unknown, whereas all the other terms are taken to the right-hand-side. The superscript k is from the previous Newton iteration. The terms in black are from the mass matrix, in blue are the nonlinear terms, in red the diffusive terms, in green the remaining terms from the stiffness matrix and in yellow the source terms.

Susceptible (Equation 15):

Ωe ∇N a n k ν s ∇N b dΩ (61) 

Epidemiologia matemática: Estudos dos efeitos da vacinação em doenças de transmissão direta

A dengue model with a dynamic aedes albopictus vector population

Calibration of a seir-sei epidemic model to describe the zika virus outbreak in brazil

Global analysis of an hiv/aids epidemic model

Modeling the sars epidemic

Remote sensing-based time series models for malaria early warning in the highlands of ethiopia

Statistical inference in a stochastic epidemic seir model with control intervention: Ebola as a case study

Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to covid-19, mathematical analysis, and numerical study

Simulating the spread of covid-19 via spatially-resolved susceptible-exposed-infected-recovered-deceased (seird) model with heterogeneous diffusion

A simple model for covid-19

Modelling the covid-19 epidemic and implementation of population-wide interventions in italy

A simulation of a covid-19 epidemic based on a deterministic seir model

Spreading of covid-19 in brazil: Impacts and uncertainties in social distancing strategies. medRxiv

Mathematical models in epidemiology

Mathematical models in population biology and epidemiology

Numerical simulation of a susceptible-exposedinfectious space-continuous model for the spread of rabies in raccoons across a realistic landscape

Compartmental modelling in chemical engineering: A critical review

Partial differential equations in ecology: spatial interactions and population dynamics

Modeling infectious diseases in humans and animals

Spread and dynamics of the covid-19 epidemic in italy: Effects of emergency containment measures

A multi-scale model of virus pandemic: Heterogeneous interactive entities in a globally connected world

Galerkin methods for a model of population dynamics with nonlinear diffusion

Spatial ecology via reaction-diffusion equations

Computational Grids, Generation, Adaptation and Solution Strategies

Modeling future spread of infections via mobile geolocation data and population dynamics. an application to covid-19 in brazil

The effect of human mobility and control measures on the covid-19 epidemic in china

deal. ii-a general-purpose object-oriented finite element library

The fenics project version 1.5. Archive of Numerical Software

Grins: a multiphysics framework based on the libmesh finite element library

Moose: A parallel computational framework for coupled systems of nonlinear equations

Residual-based variational multiscale 2d simulation of sediment transport with morphological changes

A new convected level-set method for gas bubble dynamics

EPIDEMIC -Epidemiology Educational Code

On the definition and the computation of the basic reproduction ratio r0 in models for infectious diseases in heterogeneous populations

libmesh: a c++ library for parallel adaptive mesh refinement/coarsening simulations

In situ visualization and data analysis for turbidity currents simulation

Dfanalyzer: Runtime dataflow analysis tool for computational science and engineering applications

A posteriori error estimation in finite element analysis

A posteriori error analysis and adaptive processes in the finite element method: Part i -error analysis

Reduced Order Methods for Modeling and Computational Reduction

Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control

Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations

This research was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001 and CAPES TecnoDigital Project 223038.014313/2020-19. This research has also received funding from CNPq and FAPERJ. We are indebted to Prof. Americo Cunha Jr., Prof. Regina Almeida and Prof. Sandra Malta for fruitful discussions and invaluable help in the understanding of epidemiological models.

A Implementation of the generic spatio-temporal SEIRD modelWe implement the generic SEIRD model similar to the EPIDEMIC software. We have used the BDF2 time discretization method, Newton's method for the nonlinear terms, and we simplify the number of the living population by considering the previous time-step solution. For all test cases the nonlinear tolerance for Newton's method is set to 10 −8 and the linear solver tolerance is set to 10 −10 . The linear solver is GMRES with ILU(0) preconditioner.