Biblio http://www.pe.tamu.edu/blasingame/data/z_zCourse_Archive/P620_08C/P620_08C_Lec_(for_class)/P620_Mod3_FunFld_07_DifEq_MlPhs_(for_class).pdf
For the Muskat's Material Balance equation:
$$\newcommand{\numD}{{\dfrac {S_o}{B_o B_g} \dfrac {dR_s}{dP} +
\dfrac {S_o}{B_o} \dfrac {k_g}{k_o} \dfrac {\mu_o}{\mu_g} \dfrac {dB_o}{dP} +
(1 - S_o - S_w) \dfrac {1}{B_g} \dfrac {dB_g}{dP} }}$$
$$\dfrac {dS_o}{dP} = \dfrac {\numD} {1 + \dfrac {k_g}{k_o} \dfrac {\mu_o}{\mu_g} }$$
All the terms on the right side are function of pressure ($P$) and saturation ($S$), the equation could be reduced to:
$$ \dfrac{dS}{dP} = f(P, S)$$
A first-order ordinary diferential equation (ODE) with the general analytical solution:
$$y(x) = C_1e^x-2x-2$$
Given the following ODE that represents a particular reservoir: $$ \dfrac{dS}{dP} = (2P + S)$$
We wish to find the first saturation values if the initial conditions are:
$P_o$ = 0, and $S_o$ = 1.
The analytical solution for the differential equation, given the initial conditions, is: $$S(P) = 3e^P-2P-2$$
Since we know the analytical solution, we will add a method getExactSolution
to return the exact values at any given P. Then we will compare the exact solution versus the numerical solutions from the available ODE solvers in rODE
.
# the ODE object library(rODE) library(ggplot2) # class declaration setClass("MuskatODE", slots = c( stack = "environment" # environment object inside the class ), contains = c("ODE") ) # Initialization method setMethod("initialize", "MuskatODE", function(.Object, ...) { .Object@stack$n <- 0 return(.Object) }) # The exact solution method setMethod("getExactSolution", "MuskatODE", function(object, t, ...) { # analytical solution return(3 * exp(t) - 2 *t - 2) # constant C1 = 3 }) # obtain the state of the ODE setMethod("getState", "MuskatODE", function(object, ...) { object@state # return the state }) # the differential equation is entered here. setMethod("getRate", "MuskatODE", function(object, state, ...) { object@rate[1] <- state[1] + 2 * state[2] # 2P + S object@rate[2] <- 1 # dP/dP object@stack$n <- object@stack$n + 1 # add 1 to the rate count object@rate # return the rate }) # constructor MuskatODE <- function(P, S) { .MuskatEuler <- new("MuskatODE") .MuskatEuler@state[1] = S # S state[1] is saturation .MuskatEuler@state[2] = P # P = t state[2] is pressure return(.MuskatEuler) }
After we define the differential equation in shape of code, or ODE object, we move to build the iterative routine for the differential equation. We start with the given initial conditions, $P =0$, and $S=1$, stopping the iterations at $P = 0.5$.
Since we don't know the size of the infinitesimal step $\Delta P$ to use in the solver, we are given three different step sizes 0.2, 0.1, 0.05
. Depending of the accuracy of the solution, at the end a satisfactory step size will be selected.
As we go calculating the numerical solution with the Euler method, we are also getting the exact result from the analytical solution $S(P) = 3e^P-2P-2$, using both to find the error and the relative error. The number of steps that the ODE solver takes is also recorded.
# application that uses the Muskat ODE solver above MuskatEulerApp <- function(stepSize) { ode <- MuskatODE(0, 1) # initial state S(0) = 1; P = 0 ode_solver <- Euler(ode) # use the Euler ODE solver ode_solver <- setStepSize(ode_solver, stepSize) rowVector <- vector("list") # calculation will be added to rowVector pres <- 0 # P = 0 i <- 1 # index of the iterations while (pres < 0.5) { state <- getState(ode_solver@ode) pres <- state[2] error <- getExactSolution(ode_solver@ode, pres) - state[1] rowVector[[i]] <- list(step_size = stepSize, # vector with calculations P = pres, S = state[1], exact = getExactSolution(ode_solver@ode, pres), error = error, rel_err = error / getExactSolution(ode_solver@ode, pres), steps = ode_solver@ode@stack$n ) ode_solver <- step(ode_solver) # advance solver one step i <- i + 1 # advance iteration } data.table::rbindlist(rowVector) # results in data table }
After the iteration model is completed, it is time to analyze the results stored in the table. We do this by sending the three step sizes for $\Delta P$ using the function lapply
. Since we don't want to see all the intermediate steps, we will limit the result table to display only the solutions at $P$ equal 0.10, 0.20, 0.30, 0.40, 0.5
.
The table has these columns or variables:
step_size
: size of the step for the solver.P
: pressureS
: saturationexact
: exact value from the analytical solutionerror
: difference between the analytical solution and the numerical solutionrel_err
: relative errorsteps
: number of steps taken by the solver# get a summary table for different step sizes get_table <- function(stepSize) { dt <- MuskatEulerApp(stepSize) dt[round(P, 2) %in% c(0.10, 0.20, 0.30, 0.40, 0.5)] } # vector with some step sizes step_sizes <- c(0.2, 0.1, 0.05) dt_li <- lapply(step_sizes, get_table) # create a data table data.table::rbindlist(dt_li) # bind the data tables
We can see from the table that for the smaller step size of 0.05
at $P = 0.5$, we get a relative error of approximately 0.0306
, or 3%. We willimprove this by using improved solvers.
MuskatRK4App <- function(stepSize) { ode <- MuskatODE(0, 1) ode_solver <- RK4(ode) ode_solver <- setStepSize(ode_solver, stepSize) rowVector <- vector("list") pres <- 0 i <- 1 while (pres < 0.5) { state <- getState(ode_solver@ode) pres <- state[2] error <- getExactSolution(ode_solver@ode, pres) - state[1] rowVector[[i]] <- list(step_size = stepSize, P = pres, S = state[1], exact = getExactSolution(ode_solver@ode, pres), error = error, rel_err = error / getExactSolution(ode_solver@ode, pres), steps = ode_solver@ode@stack$n ) ode_solver <- step(ode_solver) i <- i + 1 } data.table::rbindlist(rowVector) } # get a summary table for different step sizes get_table <- function(stepSize) { dt <- MuskatRK4App(stepSize) dt[round(P, 2) %in% c(0.10, 0.20, 0.30, 0.40, 0.5)] } step_sizes <- c(0.2, 0.1, 0.05) dt_li <- lapply(step_sizes, get_table) data.table::rbindlist(dt_li)
We see above that we are repeating code when the only parameter needed to be changed is the ODE solver. In cases where we want to test different ODE solvers it is more convenient to use the function ODESolverFactory
and send the solver (Euler, RK4, Verlet, etc.), as a parameter.
ComparisonMuskatODEApp <- function(solver, stepSize) { ode <- MuskatODE(0, 1) solver_factory <- ODESolverFactory() ode_solver <- createODESolver(solver_factory, ode, solver) ode_solver <- setStepSize(ode_solver, stepSize) rowVector <- vector("list") pres <- 0 i <- 1 while (pres < 0.5001) { state <- getState(ode_solver@ode) pres <- state[2] error <- getExactSolution(ode_solver@ode, pres) - state[1] rowVector[[i]] <- list(solver = solver, step_size = stepSize, P = pres, S = state[1], exact = getExactSolution(ode_solver@ode, pres), error = error, rel_err = error / getExactSolution(ode_solver@ode, pres), steps = ode_solver@ode@stack$n ) ode_solver <- step(ode_solver) pres <- pres + getStepSize(ode_solver) # step size retrievd from ODE solver i <- i + 1 } data.table::rbindlist(rowVector) } # get a summary table for different step sizes create_table <- function(stepSize, solver) { dt <- ComparisonMuskatODEApp(solver, stepSize) if (!solver == "RK45") dt[round(P, 2) %in% c(0.10, 0.20, 0.30, 0.40, 0.5)] else dt }
# Create summary table for ODE solver Euler step_sizes <- c(0.2, 0.1, 0.05) dt_li <- lapply(step_sizes, create_table, solver = "Euler") data.table::rbindlist(dt_li)
# Create summary table for ODE solver EulerRichardson step_sizes <- c(0.2, 0.1, 0.05) dt_li <- lapply(step_sizes, create_table, solver = "EulerRichardson") data.table::rbindlist(dt_li)
# Create summary table for ODE solver RK4 step_sizes <- c(0.2, 0.1, 0.05) dt_li <- lapply(step_sizes, create_table, solver = "RK4") data.table::rbindlist(dt_li)
# Create summary table for ODE solver RK45 step_sizes <- c(0.2, 0.1, 0.05) # do not round because RK45 makes variable step sizes dt_li <- lapply(step_sizes, create_table, solver = "RK45") data.table::rbindlist(dt_li)
And so on. We could this better. We will show how to do this using nested lapply
functions.
# vectors for the solvers and step sizes step_sizes <- c(0.2, 0.1, 0.05) solvers <- c("Euler", "EulerRichardson", "Verlet", "RK4", "RK45")
# nested lapply to iterate through solvers and step sizes df_li <- lapply(solvers, function(svr) lapply(step_sizes, function(stepsz) create_table(stepsz, svr))) # join the resulting dataframes df_all <- data.table::rbindlist(unlist(df_li, recursive = FALSE)) df_all
ggplot(df_all, aes(x = P, y = rel_err, group = step_size, fill = solver )) + geom_line() + geom_area(stat = "identity") + facet_grid(step_size ~ solver)
In this last plot we want to compare the relative error of the ODE solvers that show the least error: RK4 and RK45. At the same time, we will exclude the step size of 0.2 since its error magnitude would hide those of the smaller steps.
ggplot(subset(df_all, solver %in% c("RK4", "RK45") & step_size %in% c(0.1, 0.05)), aes(x = P, y = rel_err, group = step_size, fill = solver )) + geom_line() + geom_area(stat = "identity") + facet_grid(step_size ~ solver)
It turns out that RK45
is the more accurate of all the solvers tested here. The relative error at a step value of 0.05 is in the order of 1E-8 to 1 E-12, while RK4 for the same step size, the relative error ranges between 1E-8 to 9E-8; still good. In fact, we could choose with some peace of mind RK4, unless the response is unstable that merits the switch to an adaptive step solver such as RK45. RK4 is widely used by its balance of accuracy and computational speed.
Since the conditions in the equation have not change, we will use the same ODE object above MuskatODE
.
$$ \dfrac{dS}{dP} = (2P + S)$$
We are given the pressure value of $P = 3$ and the step size of 0.05. We would need to give only minor touches to the ODE solver algorithm.
This time, we will add an additional parameter to the Muskat application: the pressure of interest or pmax
, so could solve similar problems.
MuskatODEApp <- function(solver, stepSize, pmax) { ode <- MuskatODE(0, 1) solver_factory <- ODESolverFactory() ode_solver <- createODESolver(solver_factory, ode, solver) ode_solver <- setStepSize(ode_solver, stepSize) rowVector <- vector("list") pres <- 0 i <- 1 while (pres < pmax) { state <- getState(ode_solver@ode) pres <- state[2] error <- getExactSolution(ode_solver@ode, pres) - state[1] rowVector[[i]] <- list(solver = solver, step_size = stepSize, P = pres, S = state[1], exact = getExactSolution(ode_solver@ode, pres), error = error, rel_err = error / getExactSolution(ode_solver@ode, pres), steps = ode_solver@ode@stack$n ) ode_solver <- step(ode_solver) pres <- pres + getStepSize(ode_solver) # step size retrievd from ODE solver i <- i + 1 } data.table::rbindlist(rowVector) } solver <- "RK4" stepSize <- 0.05 pmax <- 3.0 dt <- MuskatODEApp(solver, stepSize, pmax) dt
last_row <- tail(dt, 1) last_row
ggplot(dt, aes(x = P, y = S)) + geom_point()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.