ode.py 4.69 KB
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1
import os
2
import math
Dion Haefner's avatar
Dion Haefner committed
3 4 5 6 7 8 9 10
import numpy as np
from scipy import integrate
import matplotlib as mpl
mpl.use("agg")
import matplotlib.pyplot as plt

from dorie.utilities.vtktools import pvd_reader
from dorie.testtools.utilities.decorators import evaluation_script
11 12
from dorie.testtools.utilities.read_yml import read_yml, get_parameter_data
from dorie.testtools.utilities.read_h5 import read_h5_dataset
Dion Haefner's avatar
Dion Haefner committed
13 14 15 16
from dorie.testtools.utilities import richards_equation

@evaluation_script
def evaluate(iniinfo,runtime):
Santiago Ospina's avatar
Santiago Ospina committed
17
    pvd = iniinfo["richards.output.fileName"] + ".pvd"
Dion Haefner's avatar
Dion Haefner committed
18 19 20
    vtks = pvd_reader.read(pvd)
    vtk_file = vtks[-1]

21
    influx = float(iniinfo["_ode.flux"])
Santiago Ospina's avatar
Santiago Ospina committed
22
    initial_head = float(iniinfo["richards.initial.headLower"])
Dion Haefner's avatar
Dion Haefner committed
23 24
    extensions = list(map(float,iniinfo["grid.extensions"].split()))

25
    # get parameter data
26
    param_data = read_yml(iniinfo["richards.parameters.file"])
27
    mapping_file = iniinfo["grid.mapping.file"]
28 29 30 31 32 33 34 35 36 37 38 39 40 41

    # homogeneous setting
    if mapping_file in ["None", "none", ""]:
        medium_index = int(iniinfo["grid.globalIndex"])
        param = get_parameter_data(medium_index, volumes=param_data["volumes"])
        if param["type"] != "MvG":
            raise NotImplementedError("This testing function only supports "
                                      "Mualem van Genuchten parameterizations. "
                                      "Parameterization was: {}".format(
                                        param["type"]))

        parfun = lambda y: param["parameters"]

    # 1D heterogeneity
42

43
    else:
44
        mapping_dataset = iniinfo["grid.mapping.volume"]
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
        mapping = read_h5_dataset(mapping_file, mapping_dataset)

        # check mapping shape
        if np.any(np.array(mapping.shape[1:]) > 1):
            raise RuntimeError("ODE solutions only work for pseudo 1D cases, "
                               "with a single cell in horizontal directions.")

        def parfun(y):
            # cell index
            cell = math.floor(y * mapping.shape[0] / extensions[-1])
            index = mapping[cell, ...].flatten()
            param = get_parameter_data(index, volumes=param_data["volumes"])
            if param["type"] != "MvG":
                raise NotImplementedError("This testing function only supports "
                                      "Mualem van Genuchten parameterizations. "
                                      "Parameterization was: {}".format(
                                        param["type"]))
            return param["parameters"]
Dion Haefner's avatar
Dion Haefner committed
63 64 65

    print("Operating on file: {}".format(vtk_file))

66
    grid = vtk_file.read(["head","flux","K_0"])
Dion Haefner's avatar
Dion Haefner committed
67
    cells = grid.cellCenters()
68 69
    y = cells[...,1].flatten()
    tri, connectivity = grid.triangulation()
Dion Haefner's avatar
Dion Haefner committed
70 71
    data = grid.data()
    iy = np.argsort(y)
72
    iy_inv = np.arange(len(y))[iy]
Dion Haefner's avatar
Dion Haefner committed
73 74
    y = y[iy]

75
    k0_dorie = data["K_0"].flatten()[iy]
Dion Haefner's avatar
Dion Haefner committed
76 77 78 79
    k0_calc = np.array([parfun(yi)["k0"] for yi in y]).flatten()
    maxdiff = np.max(np.abs(k0_dorie-k0_calc))
    if maxdiff > 1e-10:
        print("Encountered unexpected parameter values. Maximum difference: {:.2e}".format(maxdiff))
80 81 82 83
        return False
    del k0_dorie
    del k0_calc

Dion Haefner's avatar
Dion Haefner committed
84
    r = integrate.ode(richards_equation.richards)
85
    r.set_integrator("dopri5", atol=1e-10, rtol=1e-10)
Dion Haefner's avatar
Dion Haefner committed
86
    r.set_initial_value(initial_head, 0)
87
    r.set_f_params(influx,parfun)
Dion Haefner's avatar
Dion Haefner committed
88

89 90
    r_y = []
    for t in y:
91
        if t - r.t > 1E-7: # ODE solver can't handle step sizes of ~ 0
Dion Haefner's avatar
Dion Haefner committed
92
            r.integrate(t)
93
        r_y.append(r.y[0])
94
    r_y = np.array(r_y)
Dion Haefner's avatar
Dion Haefner committed
95
    head = data["head"].flatten()[iy]
96
    res_head = r_y - head
97 98
    average = lambda v: grid.integrate(v) / grid.integrate(1)
    l2_head = np.sqrt(average(np.square(res_head[iy_inv])))
99
    print("head L2 error: {:.2e}".format(float(l2_head)))
Dion Haefner's avatar
Dion Haefner committed
100

101
    plt.figure()
102 103
    plt.plot(y,r_y,label="ODE solution")
    plt.plot(y,head,label="DORiE solution")
104
    plt.legend(loc="best")
Santiago Ospina's avatar
Santiago Ospina committed
105
    plt.savefig("{}/head_solutions.png".format(iniinfo["richards.output.outputPath"]))
106 107 108 109
    plt.close()

    plt.figure()
    plt.plot(y,res_head)
Santiago Ospina's avatar
Santiago Ospina committed
110
    plt.savefig("{}/head_residual.png".format(iniinfo["richards.output.outputPath"]))
111 112 113 114
    plt.close()

    flux = data["flux"][iy,:]
    res_flux = flux - np.array([0.,influx,0.])[np.newaxis,:]
115
    l2_flux = [np.sqrt(average(np.square(r[iy_inv]))) for r in res_flux.T]
116
    print("flux L2 errors: {:.2e}, {:.2e}, {:.2e}".format(*l2_flux))
117 118

    plt.figure()
119
    plt.plot(y,res_flux)
Santiago Ospina's avatar
Santiago Ospina committed
120
    plt.savefig("{}/flux_residual.png".format(iniinfo["richards.output.outputPath"]))
121
    plt.close()
122

123 124
    head_tol = 1E-5 if not "_ode.head_abstol" in iniinfo else float(iniinfo["_ode.head_abstol"])
    flux_tol = abs(influx) * 1e-5 if not "_ode.flux_abstol" in iniinfo else float(iniinfo["_ode.flux_abstol"])
Dion Haefner's avatar
Dion Haefner committed
125

126
    return bool(l2_head < head_tol) and all([l < flux_tol for l in l2_flux])