A FULLY-DECOUPLED DISCONTINUOUS GALERKIN APPROXIMATION OF THE CAHN–HILLIARD–BRINKMAN–OHTA–KAWASAKI TUMOR GROWTH MODEL

. In this article, we consider the Cahn–Hilliard–Brinkman–Ohta–Kawasaki tumor growth system, which couples the Brinkman flow equations in the porous medium and the Cahn–Hilliard type equation with the nonlocal Ohta–Kawasaki term. We first construct a fully-decoupled discontinuous Galerkin method based on a decoupled, stabilized energy factorization approach and implicit-explicit Euler method in the time discretization


Introduction
Cancer is the most common malignant tumor and has become one of the highest morbidity and mortality in the world.The main difficulty in cancer detection and treatment lies in the complexity of the internal mechanisms during tumor growth.A better understanding of tumor growth can help detect and treat cancer earlier and more accurately, so it has very direct clinical significance.In the past few decades, many mathematical models have been developed and simulated to understand the mechanism of tumor growth.For examples, the nonlinear reaction-diffusion models [1][2][3][4], the hyperbolic-elliptic type models [5,6], state-parameter estimation model [7], multicellular tumor spheroid model [8] and so on.In the latest stage of tumor growth, several models are based on the hypothesis that different tissue components of the tumor are separated by a interfacial layer and therefore can be described by interface problems, it has turned out that interface models, treating the tumour as a collection of cells is a good strategy to describe the evolution and interactions of different species.If the existing mathematical models describing tumor growth are classified according to the types of interface equations, we find that there are mainly two types of interface models, namely the sharp interface model [5,6,[9][10][11][12][13], and the diffusive interface (phase-field) model [14][15][16][17][18][19][20][21][22][23]65], and the latter is considered in this paper.
In the phase-field model, the interface of the tumor is represented by a phase variable with smooth changes but a large gradient between the two distinct values.Since the phase-field model usually has some good properties such as the law of energy dissipation, it has been used to study tumor growth for more than two decades, see [5,6,10,12,19,[24][25][26] and reference therein.But it is worth noting that many developed models omit the actual interstitial fluid pressure [6,10,24,27] or just treat it as a fixed function (or constant) [13,26,28,29], although many physiological results suggest that tumor interstitial fluid pressure is a very important factor that can affect tumor prognosis, treatment, metastasis, and drug delivery, cf.[30][31][32][33][34][35][36].This prompts us to develop a comprehensive mathematical model to incorporate more details, especially the fluid pressure, to provide further insights into describing, understanding, and predicting tumor evolution.Moreover, most of the work related to the phase-field tumor growth model is devoted to modeling and theoretical analysis (e.g.well-posedness and regularity [11, 14-21, 23, 37, 38]), and there is relatively little work on algorithm design of the existing models.
In this article, inspired by the Cahn-Hilliard-Darcy/Cahn-Hilliard-Brinkman phase-field models developed in [16,20,37], we obtain the so-called Cahn-Hilliard-Brinkman-Ohta-Kawasaki tumor growth model, in which, we add the nonlocal Ohta-Kawasaki term in the Cahn-Hilliard equation and time-dependence term to the Brinkman flow model.In this way, we can address nonlocal effects to depict long-range interactions of cell species and the time-dependent behavior of interstitial fluid pressure (or velocity) at the stages of tumor growth.Then, we further consider the development of a fully discrete numerical scheme to solve this model.As mentioned above, most of the research on the phase-field tumor growth model is focused on the existence and uniqueness of weak solutions, and so there is no numerical algorithm that had been developed for this specific model.But if we expand our visions to consider other systems with similar structures to the model studied in this article, for instance, the Cahn-Hilliard equation (Ohta-Kawasaki term is absent) coupled with another type of flow field such as the incompressible Navier-Stokes equation, there exist many works focusing on the development of numerical algorithms, see [39][40][41][42][43][44][45] and reference therein.If we keep the Cahn-Hilliard equation, but further switch the flow field to the Darcy (or Brinkman) flow field, such as the Cahn-Hilliard-Darcy equation in the Hele-Shaw cell, or Cahn-Hilliard-Stokes-Darcy model, there are also some existing effective numerical schemes, see [46][47][48][49][50][51][52].
However, most of the above-mentioned numerical algorithms have some aspects that are not suitable for this particular model and the focus of this article, for example, either the algorithm only considers the time semi-discrete version assuming continuous space [44], or the nonlinear potential considered in the system is the double-well type [48], or the obtained algorithm is nonlinear and fully coupled [47].In contrast, the focus of this paper is to construct a fully discrete scheme for the Cahn-Hilliard-Brinkman-Ohta-Kawasaki system, and it is expected that the scheme can follow the linearity, decoupling, and unconditional energy stability.We adopt the discontinuous Galerkin (DG) method for spatial discretization since it offers some particular and remarkable features such as arbitrary order accuracy, local mass conservation, ready parallelization, and adaption (see [53][54][55]), etc.Meanwhile, DG might have high performance in simulating the tumor boundaries with highly irregular properties [55].It can be seen that there exists an increasing interest in applying the DG methods for solving the phase-field related models, see [56][57][58][59].However, as far as the author knows, if the logarithmic Flory-Huggins potential is used, then for either the Cahn-Hilliard-Brinkman-Ohta-Kawasaki model studied in this article or those widely concerned models like the Cahn-Hilliard-Brinkman or Cahn-Hilliard-Navier-Stokes system, how to use the DG method for the spatial discretization to obtain the decoupled and energy stable fully discrete scheme has not been resolved successfully.
Therefore, in this paper, we achieve a fully discrete DG scheme for the Cahn-Hilliard-Brinkman-Ohta-Kawasaki system with the logarithmic Flory-Huggins potential by combining the stabilized energy factorization approach with some subtle implicit-explicit treatments for nonlinear coupling terms.The scheme is highly efficient since one can efficiently solve only a sequence of elliptic equations at each time step.We also conduct rigorous error estimates, especially for the tumor interstitial fluid pressure.Some simulations related to the tumor growth process are carried out and also validated by evaluating the performance and benefits in simulating the brain tumor growth for patients based on MR (magnetic resonance) images.
The rest of this paper is organized as follows.In the next section, we give some notations and formulate the Cahn-Hilliard-Brinkman-Ohta-Kawasaki tumor growth model.In Section 3, we construct a fully decoupled DG scheme for solving the coupled system, and we also prove the energy stability of the proposed scheme.Section 4 is devoted to the error analysis, where we rigorously establish optimal error estimates for the fully discrete scheme under some suitable regularities.In Section 5, we present some numerical results to validate the developed scheme, and simulate the brain tumor growth process.Finally, some concluding remarks are given in Section 6.

Tumor growth model
In this article, we consider a modified Cahn-Hilliard phase-field model for the tumor growth, where a "growth" term is added in chemical potential, that reads as: where  denotes the tumor cell density and u is the tumor interstitial fluid velocity,  > 0 is a mobility constant, ∆ −1 stands for the inverse Laplacian operator, and φ0 = |Ω| −1 ∫︀ Ω  0 dx is the initial mass average over the domain Ω.The term ∆ −1 (︀  − φ0 )︀ with  > 0 had been used in the nonlocal Ohta-Kawasaki model to describe the phase change of diblock copolymers, see [60].In the context of tumor growth, this term works as a growth (tumor proliferation or death) term and  can be positive or negative.Namely, when  ≤ 0, the process of phase separation is enhanced, which is used to describe the enhanced proliferation of tumor cells; when  > 0, the term can suppress both the coarsening process and the phase separation, which implies the death of tumor cells.
The interstitial fluid flow in tumors also plays an important role in describing tumor growth, metastasis and treatment (see [31,36]).The transport velocity in the interstitium of tumors is the main mechanism for nutrients supply and waste removal during tumor growth and tumor cell metastasis to distant organs.In this article, we consider that the tumor interstitial fluid pressure  and velocity u are obtained by using a time-dependent Brinkman flow [61] through the porous medium, that reads as where  > 0 is the fluid viscosity,  > 0 is the Darcy drag parameter, (u) = (︁ ∇u + (∇u)

𝑇
)︁ /2 is the velocity deformation tensor.The term ∇ denotes the surface tension force with  > 0.
The physically most relevant choice of  is given by the logarithmic Flory-Huggins energy potential (see [62,63] and the references therein), namely, where  > 2 is the energy parameter, the choice of this case is necessary to ensure that it has two local minimums that allows the co-existence of two distinct phases.The system (2.3) presents a coupling system consisting of a Cahn-Hilliard-Ohta-Kawasaki equation and a time-dependent Brinkman flow equation, where the two coupling terms include the surface tension force ∇ and the advection term ∇ • (u).
The system (2.3) is usually posed by the following initial and boundary conditions: where n is the unit outward normal vector to the boundary Ω.
We fixed some notations here.For 1 ≤  ≤ ∞, let   (Ω) and  , (Ω) be the standard Lebesgue and Sobolev spaces endowed with their usual norms, respectively.Specifically, we denote by (•, •) Ω * the standard  2 -inner product equipped with the norm as where Ω * denotes the proper space as needed below.When Ω * = Ω, we also denote (•, •) Ω by (•, •) in short.We also introduce a Hilbert space  −1 (Ω), which is the dual of the Sobolev space  1 (Ω).Recall that, if  ∈  2 (Ω) and ∫︀  = 0, there is a unique  such that We introduce the notation  = (−∆) −1 .Since (−∆) −1 is a positive self-adjoint operator and fractional power of it is well defined.The space is the completion of the space of smooth functions in the norm . The inner product is given by for ,  belonging to  −1 (Ω).If  ∈  −1 (Ω) and  ∈  2 (Ω), then we have We consider the of weak formulation of the Cahn-Hilliard-Brinkman-Ohta-Kawasaki system (2.3), which reads as follows: find (, , , u, ) such that for all (, , , v, ) (2.7) Note that the coupled system (2.3) satisfies the law of energy dissipation.To see that, we consider the following energy functional It is straightforward to show that the system (2.7) admits the law of energy, we state the result as a lemma here.
(2.9) Thus, taking the summation of five equations in (2.9), we can easily get (2.10) where we use the fact that and .
After dealing with the above equation (2.10), the desired result is obtained.The proof of Lemma 2.1 is finished.

A fully-decoupled DG scheme
Let  be a positive integer and 0 =  0 <  1 < • • • <   =  be an uniform partition of time interval [0,  ] with time step  =  +1 −   ,  = 0, 1, • • • ,  − 1. Suppose ℰ ℎ = {} is a family of non-overlap subdivision of Ω parameterized by ℎ > 0, where ℎ denotes the discrete spatial mesh size and triangle  stands for physical computation element, we define ℎ  = diam() and ℎ = max ∈ℰ ℎ ℎ  .Let Γ ℎ denote the set of all interior edges of ℰ ℎ .Let  be a segment of  shared by two elements   1 and   2 which are neighbors, associated with  once and for all, a unit normal vector n  oriented from   1 to   2 .We define formally the average and jump of a scalar or vector valued function  for interior edges and by By convention we can extend the definitions of average and jump to sides that belong to boundary Ω, as follows Some discontinuous/broken Sobolev spaces (see [55]) on the decomposition ℰ ℎ shall be recalled.For some nonnegative integer  ≥ 0, we define equipped with the broken Sobolev norm: and in particular, we will use the broken gradient semi-norm as )︃ 1/2 .
We now introduce the discontinuous finite element spaces  ℎ , X ℎ and  ℎ associated with the triangulation ℰ ℎ used in this article, define where P  () denotes the set of all polynomials of degree less than or equal to (≥ 1) on the element .
The DG discretization of the operators −∆, −∇ • [2(v)] and −∆ are enforced by the bilinear forms   ,  ℐ and   as follows: Recalling that, since we use the symmetric bilinear forms, the penalty parameter   has to be chosen large enough.
In addition, the discretization of the term ∇ is done by the bilinear form   as The DG discretization of advection term ∇ • (v) and interface term −∇, respectively, are given by In a convenient manner, here we specially design our numerical scheme to satisfy ℬ  (, u, ) = ℬ ℐ (, , u) for any ,  and u.
Based on the discontinuous Sobolev spaces and the DG operators, the energy norms of relevant spaces are defined by Whenever   (, ) ≥ 0,  ℐ (, ) ≥ 0 and   (, ) ≥ 0 (see Lem.Let Ŝℎ be the dual of  ℎ , we also introduce the inverse discrete Laplace operator (−∆ ℎ ) −1 : Ŝℎ → Ŝℎ as follows (see [64]): where   is defined as above.Let ,  ∈ Ŝℎ , we define the following inner product by and induced the norm as Consequently, for all ,  ∈ Ŝℎ and  ∈  ℎ , we have Furthermore, for all  ∈ Ŝℎ , the following Poincaré-type estimate holds where  > 0 is the usual Poincaré constant.Due to the strong nonlinearity of the energy potential function, a challenging issue to solve the system (2.3) numerically is how to design efficient schemes that preserve the energy stability of the discrete system.In this study, following the work in [62], we regularize the logarithmic bulk potential by a  2 piecewise function.More precisely, for any 0 <  ≪ 1, the regularized free energy is Then, the derivative of ̂︀  () is From now on, we consider the problem formulated with the substitute ̂︀  and ̂︀  .For convenience, we will omit the ̂︀ signs for both ̂︀  and ̂︀  .Now the domain for the regularized functional  () is (−∞, +∞).Hence, we do not need to worry about that any small fluctuation near the domain boundary (0, 1) of the numerical solution can cause the overflow.Here, a energy factorization approach is adapted to deal with the regularized functional  (), which can be split-up as  () =  1 () −  2 (), we denote by  1 () and  2 (): We also denote by  1 () =  ′ 1 () and  2 () =  ′ 2 ().It is easy to see that the derivative of  () can be split accordingly as  () =  1 () −  2 (), where From the definitions of  1 and  2 in (3.7), we can easy prove that Next, we are ready to construct a fully-decoupled, first-order semi-discrete time-marching numerical scheme for solving the system (2.3).Given the initial conditions (︀  0 ,  0 ,  0 , u 0 ,  0 )︀ , having computed (  ,   ,   , u  ,   ).We update (︀  +1 ,  +1 ,  +1 , u +1 ,  +1 )︀ for  ≥ 0 from: where u  * and  *,+1 can be calculated as Step 2. Find ũ+1 such that (3.11) Step 3. (3.12) In the above scheme, the term  (︀   ,  +1 )︀ shall be used the splitting scheme as Here, the choice of u  * can be seen in [39,40] to the interested readers.Remark 3.1.By taking the divergence for the first equation in (3.12), we obtain Now, we begin to develop the fully-decoupled DG approximation of the modified Cahn-Hilliard-Brinkman-Ohta-Kawasaki system (2.3).The fully discrete version of (3.9)-(3.14)reads as: ) Step 2. Find ũ+1 ℎ ∈ X ℎ for all v ℎ ∈ X ℎ such that Step 3. Find  +1 ℎ ∈  ℎ for all  ℎ ∈  ℎ such that We assume the initial data  0 ℎ , u 0 ℎ and  0 ℎ are good approximations of  0 , u 0 and  0 , respectively.For example, we choose  0 ℎ =  ℎ  0 , u 0 ℎ =  ℎ u 0 and  0 ℎ =  ℎ  0 , where  ℎ ,  ℎ and  ℎ denote the projection operators, respectively, see [55] for the details.
In the following, we will recall several well-known results and provide some briefly proofs for the interior penalty DG forms, which shall be used to success in performing the analysis of our scheme throughout this work.Without loss of generality, we denote by  a generic constant that is independent of  and ℎ but possibly depends on the regularity of solution.For the sake of brevity, we use the abbreviation   for the inequality  ≤ , where  > 0 denotes a generic constant independent of the mesh size, as well as the definition of  .Moreover, the constant  may indicate different values at their different appearances.Lemma 3.1 (Broken Poincaré's inequality [55]).For any  ∈ [1, +∞), the embedded inequalities hold that ‖‖   (Ω) ‖‖ DG , ‖v‖   (Ω) ‖v‖ DG , for all  ∈  ℎ and v ∈ X ℎ .Lemma 3.2 (Trace inequality [55]).For any  ∈ P  () and  ∈ , there hold Lemma 3.3 (Continuity of   ,  ℐ and   , see [55]).The bilinear forms   and  ℐ defined on  ℎ and X ℎ equipped with the energy norms are continuous, respectively.Then, we have Lemma 3.4 (Coercivity of   ,  ℐ and   , see [55]).Assume that   is sufficiently large.Then there hold Lemma 3.5 (Boundedness of   ).The bilinear form   is continuous.Moreover, the following estimate holds Proof.The proof can be seen in Appendix A. Lemma 3.6 (Boundedness of ℬ  and ℬ ℐ ).The trilinear forms ℬ  and ℬ ℐ satisfy the following results and ) It should be noted that, the operator ℬ ℐ has the same arguments.
Proof.The proof can be found in Appendix B.
Analogously to the total energy (2.8) at the continuous level, we define the discrete energy as follows The next statement, the discrete energy dissipation law, stems from a direct result of stabilized energy factorization represented in the scheme.
, with the other variables regarded as auxiliary.Then, for any ℎ and , the property of unconditional energy stability holds: 19), and using the following identity we derive 1 2 )︁ = 0. ( Based on the result (3.8), using the Taylor expansion and the facts that  ′ 1 () =  1 () and for some . Similarly, we use the the facts that )︀ , we derive that Using the definition of (3.2), we immediately have Since   is a symmetric bilinear form, using (3.1) and (3.28), we derive )︁ ≤ 0. (3.39) From (3.39) and (3.1), we have from which (3.26) follows immediately.
The discrete energy law (3.26)immediately implies the following uniform a priori estimates for   ℎ ,   ℎ , u  ℎ and   ℎ .
Then, the following estimates hold for any ℎ and : for some constant  0 > 0 that is independent of  and ℎ.

Error analysis
4.1.Error estimates for , , u and We are now in position to prove the error estimates for the tumor cell density , chemical potential , interstitial fluid velocity u and pressure .To begin with, the weak formulation of the modified Cahn-Hilliard-Brinkman-Ohta-Kawasaki system (2.3) satisfies the following truncation forms: where  +1  ,  +1 u and  +1  denote the truncation errors.In order to derive the error estimate for the numerical scheme in terms of time and space discretization, we shall assume that the weak solution to the system (2.3) is regular enough.More precisely, we assume (A) : Thus, we can easily establish the following estimates for the truncation errors, provided that the exact solutions are sufficiently smooth or in the assumption (A).
Proof.Since the proof is rather standard (see [24]), we omit the details for simplicity.
For the derivation of the error estimates, we define the projection errors and the discretization errors as follows where we denote in which the definitions of  ℎ ,  ℎ and  ℎ are given as similar in Lemma 3.8.Moreover, the operator  ℎ satisfies the following (see [55]) By subtracting (4.1)-(4.6)from the corresponding decoupled schemes (3.15)-(3.22),we derive the error equations as follows ) Theorem 4.1.Under the Assumption (A), the following inequality holds, and 1 2 and (4.17) in (4.9), respectively.Noticing the definition of (−∆ ℎ ) −1 we get and 12), we have 1 2 Applying the Cauchy-Schwarz inequality, Young's inequality and Lemma 4.1, we arrive at ) ) ) A repeated application of the Cauchy-Schwarz inequality, Hölder'inequality, Young's inequality, Lemmas 3.1 and 3.8, we get ) ) ) ) ) By virtue of the Cauchy-Schwarz inequality, Young's inequality, Lemmas 3.3 and 3.8, we obtain ) ) ) From the Cauchy-Schwarz inequality, Young's inequality, Lemma 3.8, and (3.2), (3.3), we simply bound and Since   (︀  +1 u ,  ℎ )︀ = 0, similarly, using the Lemmas 3.4 and 3.8 yields ) ) For the term  22 and  23 , from (3.13), making use of the Lagrange's mean value theorem, we give the definition of    as where  * is between ( +1 ) and  +1 ℎ .Then, from (3.5), we can derive that and Using the Cauchy-Schwarz inequality, Lemma 3.8 and (4.45), we can prove that Taking  ℎ =    in (4.7), we arrive at then, by using the Cauchy-Schwarz inequality, Young's inequality, Lemmas 3.3, 3.8 and (4.46), (4.47), the term  23 can be bounded as ) ) Now we use the Lemmas 3.6 and 3.8 to derive ) ) Combining with the above estimates, we derive 1 2 Summing (4.56) up for  = 0 to ℓ ≤  , using the fact that  0  =  0 u =  0  = 0, for  ≥ 0, using the Lemmas 3.3 and 3.4, we obtain When 0 <  ≤  0 := 1 20 < 1 0 , for any 0 < ℓ ≤  , since 1 ≤ 1 1−0 ≤ 2 and from (4.57), it can be readily seen that The application of the discrete Gronwall's inequality to (4.58), dropping some positive terms, we obtain Noting that, if we consider the sign of  ≤ 0, by Poincaré-type estimate (3.3), we can check that the above results still work.The proof is similar to the Subsection 4.2 to deal with the nonlocal Ohta-Kawasaki term.Finally, we can obtain the desired result (4.13) by applying the triangle inequality.

Improved pressure estimates
Note that the convergence order of the pressure in Theorem 4.1 is not optimal.We are now ready to derive the improved  2 error estimate for the pressure.To simplify the notation, for a sequence of functions {  ℎ }  =0 , we denote by   the increment operator Applying the increment operator   to (4.7)-(4.12),we have for  ≥ 1 For the truncation errors in (4.59)-(4.64),we have the following lemma.
Lemma 4.2.Under the Assumption (A), we have the following estimate Proof.We leave it to the interested readers since the proof is rather standard [24].
Lemma 4.3.Under the Assumption (A), the following estimates hold Proof.By setting  = 0 in (4.59)-(4.64),and taking Summing up the resulted equations, using the coercivity of   and  ℐ , applying the Cauchy-Schwarz inequality, Young's inequality, Lemmas 3.1, 3.8 and 4.2, (2.2) and (3.3).similar to the derivation of (4.49), we can obtain Then, by a simple calculation and after dropping some positive terms in (4.65), we can derive that Therefore, the first estimate in Lemma 4.3 can be easily obtained by the triangle inequality.Similarly as before, it follows from (4.57) and (4.65), one can obtain the second inequality in Lemma 4.3.
From (3.10) and (3.18), we find and and (4.82) The application of the Cauchy-Schwarz inequality, Young's inequality, Lemmas 3.5, 3.8 and Theorem 4.1 yields where we use the fact that Using (4.78) and (4.79), making use of Lemmas 3.5, 3.8 and Theorem 4.1, one obtains and

.90)
Applying the discrete Gronwall's inequality for (4.90), we get Therefore, the conclusion is completed after using the triangle inequality.
Finally, thanks to the Theorem 4.1, Lemmas 4.1-4.4,the following result can therefore be proved.
Theorem 4.2.Under the Assumption (A), we have the improved error estimate for the pressure  as  Taking the summation of (4.10) and (4.12), by means of the Cauchy-Schwarz inequality, Lemmas 3.1 and 3.5, we estimate Therefore, using the discrete inf-sup condition in Lemma 3.7, Theorem 4.1, Lemmas 4.1-4.4,(4.91) and (4.92), we get The proof is completed.

Numerical simulations
In this section, we conduct various numerical experiments to verify the theoretical results of the previous sections.Numerical simulations include tests for convergence and energy stability, as well as two patient-specific simulations of brain tumor growth, which may be of great interest for neurosurgeons because it can help them better assess the risks and benefits of surgery.

Convergence and stability tests
In the first example, we verify the accuracy and stability of the proposed algorithm.The 2D rectangular domain is set as Ω = [0, 1] × [0, 1], and the initial conditions read as: )︁ to calculate the spatial convergence rate of the proposed scheme, where the error function is defined as  ℎ  =  ℎ (x,  ) −  ℎ 2 (x,  ), as well as a similar defined for the temporal convergence rate.To verify the convergence rates of spatial errors, we take the time step size  = 0.0001 and the value  = 2, choosing the decreasing mesh size ℎ = 1/16, 1/32, 1/64, 1/128, respectively.Table 1 shows the errors and convergence rates of the velocity, phase-field variable, and pressure at  = 1, where we vary different parameter .We can see that the obtained spatial convergence rate is  (︀ ℎ 2 )︀ , which is consistent with our theoretical prediction.
To obtain the convergence rate in time, we fix the mesh size ℎ = 1/500 and choose the decreasing mesh size  = 1/10, 1/20, 1/40, 1/80.In Table 2, we show the  2 errors of the velocity, phase-field variable and pressure at  = 1 with different .It can be seen that the convergence rate is () for all variables, which is consistent with the theoretical predictions given in the previous section.
Finally, we perform energy stability tests using different time steps.In Figure 1, we plot six energy evolution curves calculated by using the time steps ranging from  = 0.5 to  = 0.001 and mesh size ℎ = 1/100, where we take  = −5 (The choice of other value  gives the similar monotonic results, we omit it here).All the obtained curves display very good monotonic attenuation, which are in accordance with the discrete energy dissipation law.

Simulations of tumor growth
In this example, we provide some numerical simulations to show the effective of nonlocal terms on the tumor growth.We consider the rectangular domain Ω = [0, 2] × [0, 2], and the initial condition for the phase , is taken to be   Figure 2 show the numerical solutions of tumor phase field  at different time levels with  = 5,  = 0,  = −5 and  = −10, respectively.We can see that, when  > 0, the tumor grows very slowly, but in some cases  ≤ 0, the tumor can grow rapidly and the patterns of tumor exhibit different shapes.In short, the introduction of nonlocal Ohta-Kawasaki term in our model can affect the tumor shapes during the tumor growth.Noting that, the choice of these parameters shows that the simulation results tally with the practical scenarios very well.

Simulations of brain tumor growth based on MR images
We mainly focus on two clinical cases of patients with brain tumors, the geometric domains and computational meshes are shown in Figure 3.The MR images of a 68-year-old male patient at three different time stages are shown in Figure 4A and 5B-5C, where we can see that the volume of the tumor increases over time significantly.The first snapshot in Figure 4A on the date of 02/2015 is chosen as the initial condition for the tumor cell density, the profile of  0 is shown in Figure 4A1, where  0 is given by  0 = 0.7 +  boundary region of the tumor, the interstitial fluid velocity is relatively high.The situation of the pressure is the opposite.We observe that the tumor interstitial fluid pressure is significantly higher than the surrounding tissues, and it is the highest at the core of the tumor.These results are consistent with the experiments given in [25,33].This is because the quickly proliferating cancer cells are much more compact than the norm cells, which pushes healthy tissue outward to form a boundary to trap interstitial fluid and pressure inside the tumor.The tumor cells always migrate or proliferate towards directions of lower interstitial pressure, and as a result, the dividing cells inside the tumor successively push other cells towards the peripheral region with the highest pressure gradient.7C, which also illustrates the effectiveness of the proposed model and numerical scheme.

Concluding remarks
In this paper, we consider the Cahn-Hilliard-Brinkman-Ohta-Kawasaki system with the Flory-Huggins logarithmic potential for brain tumor growth, and design a fully-decoupled DG method for this hydrodynamically  coupled system.The time discretization is based on a stabilized energy factorization approach and some subtle explicit-implicit treatments for nonlinear coupling terms, which can efficiently preserve the energy stability of the system.The optimal error estimates for the proposed scheme are also derived.To confirm the efficiency and accuracy of the proposed scheme, we have performed various numerical experiments and the simulation results are in good agreement with theoretical prediction and real MR images of patients.For future work, we will consider the environmental factors including nutrients and oxygen concentration in the modeling work to improve the simulation and prediction of tumor growth in the quantitative way.

Lemma 4 . 1 .
Under the Assumption (A), the truncation errors    ,   u and    satisfy

Figure 1 .
Figure 1.The time evolution of the discrete energy  ℎ with different time steps.
In the following numerical examples, we carry out the simulations of brain tumor growth.The geometric domain and meshes for the patient-specific simulations of brain tumor are constructed based on segmentation of MR images.The final time of brain tumor growth is hard to handle, so we only simulate the process of tumor growth compared with the MRI results.We set the model parameters as  = −5,  = 0.1, ℎ = 1/100,  = 0.01,  = 0.1,  = 0.02,  = 0.1.

Figure 2 .
Figure 2. The numerical solutions of tumor phase field  at different time levels with different .The graphs are arranged column-wise.

Figure 3 .
Figure 3.The geometric domains and computational meshes for the two clinical cases of patients, respectively.

Figure 4 .
Figure 4.The initial condition for the tumor cell density.(A) The MR image on the date of 02/2015; and (A1) The profile of  0 .

Figure 5 .
Figure 5.The first column shows the real MR images, the dates are on: (B) 06/2015; (C) 09/2015.The second column lists the realistic simulation results of brain tumor growth (B1-C1).The third column (B2-C2) shows the tumor interstitial fluid velocity.The last column (B3-C3) shows the interstitial fluid pressure.The graphs are arranged column-wise.

Figure 6 .
Figure 6.The initial condition for the tumor cell density.(A) The MR image on the date of 08/2018; and (A1) The profile of  0 .

Table 1 .
Numerical errors and convergence rates for (, u, ) in spatial direction with different  at  = 1.

Table 2 .
Numerical errors and convergence rates for (, u, ) in temporal direction with different  at  = 1.