diff --git a/oggm/core/flowline.py b/oggm/core/flowline.py index 54e4b3d6a..9f955bdde 100644 --- a/oggm/core/flowline.py +++ b/oggm/core/flowline.py @@ -16,6 +16,7 @@ import numpy as np import shapely.geometry as shpg import xarray as xr +from scipy.linalg import solve_banded # Optional libs try: @@ -2093,6 +2094,322 @@ def step(self, dt): return dt +class SemiImplicitModel(FlowlineModel): + """Semi implicit flowline model. + + It solves the same equation as the FluxBasedModel, but the ice flux q is + implemented as q^t = D^t * (ds/dx)^(t+1). + + It supports only a single flowline (no tributaries) with bed shapes + rectangular, trapezoidal or a mixture of both. + """ + + def __init__(self, flowlines, mb_model=None, y0=0., glen_a=None, fs=0., + inplace=False, fixed_dt=None, cfl_number=0.6, min_dt=None, + **kwargs): + """Instantiate the model. + + Parameters + ---------- + flowlines : list + the glacier flowlines + mb_model : MassBalanceModel + the mass balance model + y0 : int + initial year of the simulation + glen_a : float + Glen's creep parameter + fs : float + Oerlemans sliding parameter + inplace : bool + whether or not to make a copy of the flowline objects for the run + setting to True implies that your objects will be modified at run + time by the model (can help to spare memory) + fixed_dt : float + set to a value (in seconds) to prevent adaptive time-stepping. + cfl_number : float + For adaptive time stepping (the default), dt is chosen from the + CFL criterion (dt = cfl_number * dx^2 / max(D/w)). + Can be set to higher values compared to the FluxBasedModel. + Default is 0.6, but need further investigation. + min_dt : float + Defaults to cfg.PARAMS['cfl_min_dt']. + At high velocities, time steps can become very small and your + model might run very slowly. In production, it might be useful to + set a limit below which the model will just error. + kwargs : dict + Further keyword arguments for FlowlineModel + + """ + + super(SemiImplicitModel, self).__init__(flowlines, mb_model=mb_model, + y0=y0, glen_a=glen_a, fs=fs, + inplace=inplace, **kwargs) + + if len(self.fls) > 1: + raise ValueError('Implicit model does not work with ' + 'tributaries.') + + # convert pure RectangularBedFlowline to TrapezoidalBedFlowline with + # lambda = 0 + if isinstance(self.fls[0], RectangularBedFlowline): + self.fls[0] = TrapezoidalBedFlowline( + line=self.fls[-1].line, dx=self.fls[-1].dx, + map_dx=self.fls[-1].map_dx, surface_h=self.fls[-1].surface_h, + bed_h=self.fls[-1].bed_h, widths=self.fls[-1].widths, + lambdas=0, rgi_id=self.fls[-1].rgi_id, + water_level=self.fls[-1].water_level, gdir=None) + + if isinstance(self.fls[0], MixedBedFlowline): + if ~np.all(self.fls[0].is_trapezoid): + raise ValueError('Implicit model only works with a pure ' + 'trapezoidal flowline! But different lambdas ' + 'along the flowline possible (lambda=0 is' + 'rectangular).') + elif not isinstance(self.fls[0], TrapezoidalBedFlowline): + raise ValueError('Implicit model only works with a pure ' + 'trapezoidal flowline! But different lambdas ' + 'along the flowline possible (lambda=0 is' + 'rectangular).') + + if cfg.PARAMS['use_kcalving_for_run']: + raise NotImplementedError("Calving is not implemented in the" + "SemiImplicitModel! Set " + "cfg.PARAMS['use_kcalving_for_run'] = " + "False.") + + self.fixed_dt = fixed_dt + if min_dt is None: + min_dt = cfg.PARAMS['cfl_min_dt'] + self.min_dt = min_dt + + if cfl_number is None: + cfl_number = cfg.PARAMS['cfl_number'] + if cfl_number < 0.1: + raise InvalidParamsError("For the SemiImplicitModel you can use " + "cfl numbers in the order of 0.1 - 0.6 " + f"(you set {cfl_number}).") + self.cfl_number = cfl_number + + # Special output + self._surf_vel_fac = (self.glen_n + 2) / (self.glen_n + 1) + + # optim + nx = self.fls[-1].nx + bed_h_exp = np.concatenate(([self.fls[-1].bed_h[0]], + self.fls[-1].bed_h, + [self.fls[-1].bed_h[-1]])) + self.dbed_h_exp_dx = ((bed_h_exp[1:] - bed_h_exp[:-1]) / + self.fls[0].dx_meter) + self.D_stag = [np.zeros(nx + 1)] + self.Amat_banded = np.zeros((3, nx)) + w0 = self.fls[0]._w0_m + self.w0_stag = (w0[0:-1] + w0[1:]) / 2 + self.rhog = (self.rho * G) ** self.glen_n + + # variables needed for the calculation of some diagnostics, this + # calculations are done with @property, because they are not computed + # on the fly during the dynamic run as in FluxBasedModel + self._u_stag = [np.zeros(nx + 1)] + self._flux_stag = [np.zeros(nx + 1)] + self._slope_stag = [np.zeros(nx + 1)] + self._thick_stag = [np.zeros(nx + 1)] + self._section_stag = [np.zeros(nx + 1)] + + @property + def slope_stag(self): + slope_stag = self._slope_stag[0] + + surface_h = self.fls[0].surface_h + dx = self.fls[0].dx_meter + + slope_stag[0] = 0 + slope_stag[1:-1] = (surface_h[0:-1] - surface_h[1:]) / dx + slope_stag[-1] = slope_stag[-2] + + return [slope_stag] + + @property + def thick_stag(self): + thick_stag = self._thick_stag[0] + + thick = self.fls[0].thick + + thick_stag[1:-1] = (thick[0:-1] + thick[1:]) / 2. + thick_stag[[0, -1]] = thick[[0, -1]] + + return [thick_stag] + + @property + def section_stag(self): + section_stag = self._section_stag[0] + + section = self.fls[0].section + + section_stag[1:-1] = (section[0:-1] + section[1:]) / 2. + section_stag[[0, -1]] = section[[0, -1]] + + return [section_stag] + + @property + def u_stag(self): + u_stag = self._u_stag[0] + + slope_stag = self.slope_stag[0] + thick_stag = self.thick_stag[0] + N = self.glen_n + rhog = self.rhog + + rhogh = rhog * slope_stag ** N + + u_stag[:] = ((thick_stag**(N+1)) * self._fd * rhogh + + (thick_stag**(N-1)) * self.fs * rhogh) + + return [u_stag] + + @property + def flux_stag(self): + flux_stag = self._flux_stag[0] + + section_stag = self.section_stag[0] + u_stag = self.u_stag[0] + + flux_stag[:] = u_stag * section_stag + + return [flux_stag] + + def step(self, dt): + """Advance one step.""" + + # Just a check to avoid useless computations + if dt <= 0: + raise InvalidParamsError('dt needs to be strictly positive') + + # read out variables from current flowline + fl = self.fls[0] + dx = fl.dx_meter + width = fl.widths_m + thick = fl.thick + surface_h = fl.surface_h + + # some variables needed later + N = self.glen_n + rhog = self.rhog + + # calculate staggered variables + width_stag = (width[0:-1] + width[1:]) / 2 + w0_stag = self.w0_stag + thick_stag = (thick[0:-1] + thick[1:]) / 2. + dsdx_stag = (surface_h[1:] - surface_h[0:-1]) / dx + + # calculate diffusivity + # boundary condition D_stag_0 = D_stag_end = 0 + D_stag = self.D_stag[0] + D_stag[1:-1] = ((self._fd * thick_stag ** (N + 2) + + self.fs * thick_stag ** N) * rhog * + (w0_stag + width_stag) / 2 * + np.abs(dsdx_stag) ** (N - 1)) + + # insert stability criterion dt <= dx^2 / max(D/w) * cfl_number + if not self.fixed_dt: + divisor = np.max(np.abs(D_stag[1:-1] / width_stag)) + if divisor > cfg.FLOAT_EPS: + cfl_dt = self.cfl_number * dx ** 2 / divisor + else: + cfl_dt = dt + + if cfl_dt < dt: + dt = cfl_dt + if cfl_dt < self.min_dt: + raise RuntimeError( + 'CFL error: required time step smaller ' + 'than the minimum allowed: ' + '{:.1f}s vs {:.1f}s. Happening at ' + 'simulation year {:.1f}, fl_id {}, ' + 'bin_id {} and max_D {:.3f} m2 yr-1.' + ''.format(cfl_dt, self.min_dt, self.yr, 0, + np.argmax(np.abs(D_stag)), + divisor * cfg.SEC_IN_YEAR)) + + # calculate diagonals of Amat + d0 = dt / dx ** 2 * (D_stag[:-1] + D_stag[1:]) / width + dm = - dt / dx ** 2 * D_stag[:-1] / width + dp = - dt / dx ** 2 * D_stag[1:] / width + + # construct banded form of the matrix, which is used during solving + # (see https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_banded.html) + # original matrix: + # Amat = (np.diag(dp[:-1], 1) + + # np.diag(np.ones(len(d0)) + d0) + + # np.diag(dm[1:], -1)) + self.Amat_banded[0, 1:] = dp[:-1] + self.Amat_banded[1, :] = np.ones(len(d0)) + d0 + self.Amat_banded[2, :-1] = dm[1:] + + # correction term for glacier bed (original equation is an equation for + # the surface height s, which is transformed in an equation for h, as + # s = h + b the term below comes from the '- b' + b_corr = - D_stag * self.dbed_h_exp_dx + + # prepare rhs + smb = self.get_mb(surface_h, self.yr, fl_id=0) + rhs = thick + smb * dt + dt / width * (b_corr[:-1] - b_corr[1:]) / dx + + # solve matrix and update flowline thickness + thick_new = utils.clip_min( + solve_banded((1, 1), self.Amat_banded, rhs), + 0) + fl.thick = thick_new + + # Next step + self.t += dt + + return dt + + def get_diagnostics(self, fl_id=-1): + """Obtain model diagnostics in a pandas DataFrame. + + Parameters + ---------- + fl_id : int + the index of the flowline of interest, from 0 to n_flowline-1. + Default is to take the last (main) one + + Returns + ------- + a pandas DataFrame, which index is distance along flowline (m). Units: + - surface_h, bed_h, ice_tick, section_width: m + - section_area: m2 + - slope: - + - ice_flux, tributary_flux: m3 of *ice* per second + - ice_velocity: m per second (depth-section integrated) + - surface_ice_velocity: m per second (corrected for surface - simplified) + """ + import pandas as pd + + fl = self.fls[fl_id] + nx = fl.nx + + df = pd.DataFrame(index=fl.dx_meter * np.arange(nx)) + df.index.name = 'distance_along_flowline' + df['surface_h'] = fl.surface_h + df['bed_h'] = fl.bed_h + df['ice_thick'] = fl.thick + df['section_width'] = fl.widths_m + df['section_area'] = fl.section + + # Staggered + var = self.slope_stag[fl_id] + df['slope'] = (var[1:nx+1] + var[:nx])/2 + var = self.flux_stag[fl_id] + df['ice_flux'] = (var[1:nx+1] + var[:nx])/2 + var = self.u_stag[fl_id] + df['ice_velocity'] = (var[1:nx+1] + var[:nx])/2 + df['surface_ice_velocity'] = df['ice_velocity'] * self._surf_vel_fac + + return df + + class FileModel(object): """Duck FlowlineModel which actually reads data out of a nc file.""" diff --git a/oggm/tests/conftest.py b/oggm/tests/conftest.py index bf92b7075..3a3b349ee 100644 --- a/oggm/tests/conftest.py +++ b/oggm/tests/conftest.py @@ -286,6 +286,23 @@ def hef_copy_gdir_base(request, test_dir): return init_hef(rgi_id='RGI50-11.99999') +@pytest.fixture(scope='module') +def hef_elev_gdir_base(request, test_dir): + """ Same as hef_gdir, but with elevation bands instead of centerlines. + Further, also a different RGI ID so that a unique directory is created + (as the RGI ID defines the final folder name) and to have no + collisions with hef_gdir + """ + try: + module = request.module + border = module.DOM_BORDER if module.DOM_BORDER is not None else 40 + return init_hef(border=border, flowline_type='elevation_bands', + rgi_id='RGI50-11.99998') + except AttributeError: + return init_hef(flowline_type='elevation_bands', + rgi_id='RGI50-11.99998') + + @pytest.fixture(scope='class') def hef_gdir(hef_gdir_base, class_case_dir): """ Provides a copy of the base Hintereisenferner glacier directory in @@ -302,3 +319,11 @@ def hef_copy_gdir(hef_copy_gdir_base, class_case_dir): """ return tasks.copy_to_basedir(hef_copy_gdir_base, base_dir=class_case_dir, setup='all') + + +@pytest.fixture(scope='class') +def hef_elev_gdir(hef_elev_gdir_base, class_case_dir): + """ Same as hef_gdir but with elevation bands instead of centerlines. + """ + return tasks.copy_to_basedir(hef_elev_gdir_base, base_dir=class_case_dir, + setup='all') diff --git a/oggm/tests/funcs.py b/oggm/tests/funcs.py index cc2f2ea7c..99082d842 100644 --- a/oggm/tests/funcs.py +++ b/oggm/tests/funcs.py @@ -159,6 +159,39 @@ def dummy_mixed_bed(deflambdas=3.5, map_dx=100., mixslice=None): return [fls] +def dummy_mixed_trap_rect_bed(deflambdas=2., map_dx=100., mixslice=None): + dx = 1. + nx = 200 + + surface_h = np.linspace(3000, 1000, nx) + bed_h = surface_h + shape = np.ones(nx) * np.NaN + is_trapezoid = ~np.isfinite(shape) + lambdas = shape * 0. + lambdas[is_trapezoid] = deflambdas + # use a second value for lambda in between + lambdas[15:20] = deflambdas / 2 + widths_m = bed_h * 0. + 10 + if mixslice: + lambdas[mixslice] = 0 + widths_m[mixslice] = 25 + else: + lambdas[0:10] = 0 + widths_m[0:10] = 25 + section = bed_h * 0. + + coords = np.arange(0, nx - 0.5, 1) + line = shpg.LineString(np.vstack([coords, coords * 0.]).T) + + fls = flowline.MixedBedFlowline(line=line, dx=dx, map_dx=map_dx, + surface_h=surface_h, bed_h=bed_h, + section=section, bed_shape=shape, + is_trapezoid=is_trapezoid, + lambdas=lambdas, widths_m=widths_m) + + return [fls] + + def dummy_trapezoidal_bed(hmax=3000., hmin=1000., nx=200, map_dx=100., def_lambdas=2): @@ -326,7 +359,24 @@ def get_test_dir(): return _TEST_DIR -def init_hef(reset=False, border=40, logging_level='INFO', rgi_id=None): +def init_hef(reset=False, border=40, logging_level='INFO', rgi_id=None, + flowline_type='centerlines'): + """ + + Parameters + ---------- + reset + border + logging_level + rgi_id + flowline_type : str + Select which flowline type should be used. + Options: 'centerlines' (default), 'elevation_bands' + + Returns + ------- + + """ from oggm.core import gis, inversion, climate, centerlines, flowline import geopandas as gpd @@ -373,15 +423,23 @@ def init_hef(reset=False, border=40, logging_level='INFO', rgi_id=None): return gdir gis.define_glacier_region(gdir) - execute_entity_task(gis.glacier_masks, [gdir]) - execute_entity_task(centerlines.compute_centerlines, [gdir]) - centerlines.initialize_flowlines(gdir) + if flowline_type == 'centerlines': + execute_entity_task(gis.glacier_masks, [gdir]) + execute_entity_task(centerlines.compute_centerlines, [gdir]) + centerlines.initialize_flowlines(gdir) + centerlines.catchment_area(gdir) + centerlines.catchment_intersections(gdir) + centerlines.catchment_width_geom(gdir) + centerlines.catchment_width_correction(gdir) + elif flowline_type == 'elevation_bands': + execute_entity_task(tasks.simple_glacier_masks, [gdir]) + execute_entity_task(tasks.elevation_band_flowline, [gdir]) + execute_entity_task(tasks.fixed_dx_elevation_band_flowline, [gdir]) + else: + raise NotImplementedError(f'flowline_type {flowline_type} not' + f'implemented!') centerlines.compute_downstream_line(gdir) centerlines.compute_downstream_bedshape(gdir) - centerlines.catchment_area(gdir) - centerlines.catchment_intersections(gdir) - centerlines.catchment_width_geom(gdir) - centerlines.catchment_width_correction(gdir) climate.process_custom_climate_data(gdir) if rgi_id is not None: gdir.rgi_id = 'RGI50-11.00897' @@ -419,7 +477,8 @@ def to_optimize(x): write=True) inversion.filter_inversion_output(gdir) - inversion.distribute_thickness_interp(gdir, varname_suffix='_interp') + if flowline_type == 'centerlines': + inversion.distribute_thickness_interp(gdir, varname_suffix='_interp') inversion.distribute_thickness_per_altitude(gdir, varname_suffix='_alt') flowline.init_present_time_glacier(gdir) diff --git a/oggm/tests/test_models.py b/oggm/tests/test_models.py index f80784092..0f9b38ee9 100644 --- a/oggm/tests/test_models.py +++ b/oggm/tests/test_models.py @@ -38,7 +38,7 @@ flowline_from_dataset, FileModel, run_constant_climate, run_random_climate, run_from_climate_data, equilibrium_stop_criterion, - run_with_hydro) + run_with_hydro, SemiImplicitModel) from oggm.core.dynamic_spinup import ( run_dynamic_spinup, run_dynamic_mu_star_calibration, dynamic_mu_star_run_with_dynamic_spinup, @@ -5529,3 +5529,226 @@ def test_merged_simulation(self): # Merged glacier should have a larger area after 200yrs from advancing assert (ds_entity.area.isel(time=200).sum() < ds_merged.area.isel(time=200)) + + +class TestSemiImplicitModel: + + @pytest.mark.slow + def test_equilibrium(self, hef_elev_gdir, inversion_params): + # As long as hef_gdir uses 1, we need to use 1 here as well + cfg.PARAMS['trapezoid_lambdas'] = 1 + cfg.PARAMS['downstream_line_shape'] = 'trapezoidal' + init_present_time_glacier(hef_elev_gdir) + cfg.PARAMS['min_ice_thick_for_length'] = 1 + + mb_mod = massbalance.ConstantMassBalance(hef_elev_gdir) + + fls = hef_elev_gdir.read_pickle('model_flowlines') + model = SemiImplicitModel(fls, mb_model=mb_mod, y0=0, + fs=inversion_params['inversion_fs'], + glen_a=inversion_params['inversion_glen_a'], + mb_elev_feedback='never') + + ref_vol = model.volume_km3 + ref_area = model.area_km2 + ref_len = model.fls[-1].length_m + + np.testing.assert_allclose(ref_area, hef_elev_gdir.rgi_area_km2) + + model.run_until_equilibrium(rate=1e-5) + assert model.yr >= 50 + after_vol = model.volume_km3 + after_area = model.area_km2 + after_len = model.fls[-1].length_m + + np.testing.assert_allclose(ref_vol, after_vol, rtol=0.2) + np.testing.assert_allclose(ref_area, after_area, rtol=0.01) + np.testing.assert_allclose(ref_len, after_len, atol=200.01) + + # compare to FluxBasedModel + model_flux = FluxBasedModel(fls, mb_model=mb_mod, y0=0., + fs=inversion_params['inversion_fs'], + glen_a=inversion_params['inversion_glen_a'], + mb_elev_feedback='never') + model_flux.run_until_equilibrium(rate=1e-5) + + assert model.yr == model_flux.yr + + np.testing.assert_allclose(model_flux.volume_km3, after_vol, rtol=6e-4) + np.testing.assert_allclose(model_flux.area_km2, after_area, rtol=6e-5) + np.testing.assert_allclose(model_flux.fls[-1].length_m, after_len) + + # now glacier wide with MultipleFlowlineMassBalance + cl = massbalance.ConstantMassBalance + mb_mod = massbalance.MultipleFlowlineMassBalance(hef_elev_gdir, + mb_model_class=cl) + + model = SemiImplicitModel(fls, mb_model=mb_mod, y0=0, + fs=inversion_params['inversion_fs'], + glen_a=inversion_params['inversion_glen_a'], + mb_elev_feedback='never') + model.run_until_equilibrium(rate=1e-5) + + model_flux = FluxBasedModel(fls, mb_model=mb_mod, y0=0., + fs=inversion_params['inversion_fs'], + glen_a=inversion_params['inversion_glen_a'], + mb_elev_feedback='never') + model_flux.run_until_equilibrium(rate=1e-5) + + assert model.yr >= 50 + assert model.yr == model_flux.yr + + after_vol = model.volume_km3 + after_area = model.area_km2 + after_len = model.fls[-1].length_m + + np.testing.assert_allclose(model_flux.volume_km3, after_vol, rtol=5e-4) + np.testing.assert_allclose(model_flux.area_km2, after_area, rtol=2e-5) + np.testing.assert_allclose(model_flux.fls[-1].length_m, after_len, + atol=100.1) + + @pytest.mark.slow + def test_random(self, hef_elev_gdir, inversion_params): + cfg.PARAMS['store_model_geometry'] = True + # As long as hef_gdir uses 1, we need to use 1 here as well + cfg.PARAMS['trapezoid_lambdas'] = 1 + cfg.PARAMS['downstream_line_shape'] = 'trapezoidal' + + init_present_time_glacier(hef_elev_gdir) + run_random_climate(hef_elev_gdir, nyears=100, seed=6, + fs=inversion_params['inversion_fs'], + glen_a=inversion_params['inversion_glen_a'], + bias=0, output_filesuffix='_rdn', + evolution_model=SemiImplicitModel) + run_constant_climate(hef_elev_gdir, nyears=100, + fs=inversion_params['inversion_fs'], + glen_a=inversion_params['inversion_glen_a'], + bias=0, output_filesuffix='_ct', + evolution_model=SemiImplicitModel) + + paths = [hef_elev_gdir.get_filepath('model_geometry', filesuffix='_rdn'), + hef_elev_gdir.get_filepath('model_geometry', filesuffix='_ct'), + ] + + for path in paths: + model = FileModel(path) + vol = model.volume_km3_ts() + area = model.area_km2_ts() + np.testing.assert_allclose(vol.iloc[0], np.mean(vol), + rtol=0.26) + np.testing.assert_allclose(area.iloc[0], np.mean(area), + rtol=0.1) + + @pytest.mark.slow + def test_sliding_and_compare_to_fluxbased(self, hef_elev_gdir, + inversion_params): + cfg.PARAMS['store_model_geometry'] = True + cfg.PARAMS['store_fl_diagnostics'] = True + # As long as hef_gdir uses 1, we need to use 1 here as well + cfg.PARAMS['trapezoid_lambdas'] = 1 + cfg.PARAMS['downstream_line_shape'] = 'trapezoidal' + init_present_time_glacier(hef_elev_gdir) + cfg.PARAMS['min_ice_thick_for_length'] = 1 + + start_time_impl = time.time() + run_random_climate(hef_elev_gdir, nyears=1000, seed=6, + fs=5.7e-20, + glen_a=inversion_params['inversion_glen_a'], + bias=0, output_filesuffix='_implicit_run', + evolution_model=SemiImplicitModel) + impl_time_needed = time.time() - start_time_impl + + start_time_expl = time.time() + run_random_climate(hef_elev_gdir, nyears=1000, seed=6, + fs=5.7e-20, + glen_a=inversion_params['inversion_glen_a'], + bias=0, output_filesuffix='_fluxbased_run', + evolution_model=FluxBasedModel) + flux_time_needed = time.time() - start_time_expl + + fmod_impl = FileModel( + hef_elev_gdir.get_filepath('model_geometry', + filesuffix='_implicit_run')) + fmod_flux = FileModel( + hef_elev_gdir.get_filepath('model_geometry', + filesuffix='_fluxbased_run')) + + # check that the two runs are close the whole time + assert utils.rmsd(fmod_impl.volume_km3_ts(), + fmod_flux.volume_km3_ts()) < 4e-4 + assert utils.rmsd(fmod_impl.area_km2_ts(), + fmod_flux.area_km2_ts()) < 8e-3 + + years = np.arange(0, 1001) + if do_plot: + plt.figure() + plt.plot(years, fmod_impl.volume_km3_ts(), 'r') + plt.plot(years, fmod_flux.volume_km3_ts(), 'b') + plt.title('Compare Volume') + plt.xlabel('years') + plt.ylabel('[km3]') + plt.legend(['Implicit', 'Flux'], loc=2) + + plt.figure() + plt.plot(years, fmod_impl.area_km2_ts(), 'r') + plt.plot(years, fmod_flux.area_km2_ts(), 'b') + plt.title('Compare Area') + plt.xlabel('years') + plt.ylabel('[km2]') + plt.legend(['Implicit', 'Flux'], loc=2) + + plt.show() + + for year in years: + fmod_impl.run_until(year) + fmod_flux.run_until(year) + + np.testing.assert_allclose(fmod_flux.fls[-1].length_m, + fmod_impl.fls[-1].length_m, + atol=100.1) + assert utils.rmsd(fmod_impl.fls[-1].thick, + fmod_flux.fls[-1].thick) < 1.2 + + # compare velocities + f = hef_elev_gdir.get_filepath('fl_diagnostics', + filesuffix='_implicit_run') + with xr.open_dataset(f, group='fl_0') as ds_impl: + ds_impl = ds_impl.load() + + f = hef_elev_gdir.get_filepath('fl_diagnostics', + filesuffix='_fluxbased_run') + with xr.open_dataset(f, group='fl_0') as ds_flux: + ds_flux = ds_flux.load() + + # only used to plot the velocity with the largest difference + max_velocity_rmsd = 0 + max_velocity_year = 0 + + for year in np.arange(1, 1001): + velocity_impl = ds_impl.ice_velocity_myr.loc[{'time': year}].values + velocity_flux = ds_flux.ice_velocity_myr.loc[{'time': year}].values + + velocity_rmsd = utils.rmsd(velocity_impl, velocity_flux) + if velocity_rmsd > max_velocity_rmsd: + max_velocity_rmsd = velocity_rmsd + max_velocity_year = year + + assert velocity_rmsd < 4.5 + + if do_plot: + plt.figure() + plt.plot(ds_impl.ice_velocity_myr.loc[{'time': max_velocity_year}].values, + 'r') + plt.plot(ds_flux.ice_velocity_myr.loc[{'time': max_velocity_year}].values, + 'b') + plt.title(f'Compare Velocity at year {max_velocity_year}') + plt.xlabel('gridpoints along flowline') + plt.ylabel('[m yr-1]') + plt.legend(['Implicit', 'Flux'], loc=2) + + plt.show() + + # Testing the run times should be the last test, as it is not + # independent of the computing environment and by putting it as the + # last test all other tests will still be executed + assert impl_time_needed < flux_time_needed / 2 diff --git a/oggm/tests/test_numerics.py b/oggm/tests/test_numerics.py index 525436136..c87070204 100644 --- a/oggm/tests/test_numerics.py +++ b/oggm/tests/test_numerics.py @@ -21,19 +21,21 @@ dummy_noisy_bed, dummy_parabolic_bed, dummy_trapezoidal_bed, dummy_width_bed, dummy_width_bed_tributary, bu_tidewater_bed, - dummy_bed_tributary_tail_to_head) + dummy_bed_tributary_tail_to_head, + dummy_mixed_trap_rect_bed) # after oggm.test import matplotlib.pyplot as plt from oggm.core.flowline import (KarthausModel, FluxBasedModel, MassRedistributionCurveModel, - MassConservationChecker) + MassConservationChecker, SemiImplicitModel) from oggm.tests.ext.sia_fluxlim import MUSCLSuperBeeModel FluxBasedModel = partial(FluxBasedModel, inplace=True) KarthausModel = partial(KarthausModel, inplace=True) MUSCLSuperBeeModel = partial(MUSCLSuperBeeModel, inplace=True) +SemiImplicitModel = partial(SemiImplicitModel, inplace=True) pytestmark = pytest.mark.test_env("numerics") do_plot = False @@ -60,7 +62,8 @@ def tearDown(self): @pytest.mark.slow def test_constant_bed(self): - models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel] + models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel, + SemiImplicitModel] lens = [] surface_h = [] @@ -92,41 +95,52 @@ def test_constant_bed(self): plt.plot(yrs, lens[0], 'r') plt.plot(yrs, lens[1], 'b') plt.plot(yrs, lens[2], 'g') + plt.plot(yrs, lens[3], 'm') plt.title('Compare Length') plt.xlabel('years') plt.ylabel('[m]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(yrs, volume[0], 'r') plt.plot(yrs, volume[1], 'b') plt.plot(yrs, volume[2], 'g') + plt.plot(yrs, volume[3], 'm') plt.title('Compare Volume') plt.xlabel('years') plt.ylabel('[km^3]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(fls[-1].bed_h, 'k') plt.plot(surface_h[0], 'r') plt.plot(surface_h[1], 'b') plt.plot(surface_h[2], 'g') + plt.plot(surface_h[3], 'm') plt.title('Compare Shape') plt.xlabel('[m]') plt.ylabel('Elevation [m]') - plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=3) + plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee', + 'Implicit'], loc=3) plt.show() np.testing.assert_almost_equal(lens[0][-1], lens[1][-1]) + np.testing.assert_almost_equal(lens[3][-1], lens[1][-1]) np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=3e-3) np.testing.assert_allclose(volume[1][-1], volume[2][-1], atol=3e-3) + np.testing.assert_allclose(volume[3][-1], volume[2][-1], atol=3e-3) assert utils.rmsd(lens[0], lens[2]) < 50. assert utils.rmsd(lens[1], lens[2]) < 50. + assert utils.rmsd(lens[3], lens[2]) < 50. assert utils.rmsd(volume[0], volume[2]) < 2e-3 assert utils.rmsd(volume[1], volume[2]) < 2e-3 + assert utils.rmsd(volume[3], volume[2]) < 2e-3 assert utils.rmsd(surface_h[0], surface_h[2]) < 1.0 assert utils.rmsd(surface_h[1], surface_h[2]) < 1.0 + assert utils.rmsd(surface_h[3], surface_h[2]) < 1.0 def test_length(self): @@ -187,31 +201,37 @@ def test_staggered_diagnostics(self): mb = LinearMassBalance(2600.) fls = dummy_constant_bed() - model = FluxBasedModel(fls, mb_model=mb, y0=0.) - model.run_until(700) + model_flux = FluxBasedModel(fls, mb_model=mb, y0=0.) + model_flux.run_until(700) assert_allclose(mb.get_specific_mb(fls=fls), 0, atol=10) - # Check the flux just for fun - fl = model.flux_stag[0] - assert fl[0] == 0 + model_impl = SemiImplicitModel(fls, mb_model=mb, y0=0.) + model_impl.run_until(700) - # Now check the diags - df = model.get_diagnostics() - fl = model.fls[0] - df['my_flux'] = np.cumsum(mb.get_annual_mb(fl.surface_h) * - fl.widths_m * fl.dx_meter * - cfg.SEC_IN_YEAR).clip(0) + for model in [model_flux, model_impl]: + # Check the flux just for fun + fl = model.flux_stag[0] + assert fl[0] == 0 - df = df.loc[df['ice_thick'] > 0] + # Now check the diags + df = model.get_diagnostics() + fl = model.fls[0] + df['my_flux'] = np.cumsum(mb.get_annual_mb(fl.surface_h) * + fl.widths_m * fl.dx_meter * + cfg.SEC_IN_YEAR).clip(0) - # Also convert ours - df['ice_flux'] *= cfg.SEC_IN_YEAR - df['ice_velocity'] *= cfg.SEC_IN_YEAR - df['tributary_flux'] *= cfg.SEC_IN_YEAR + df = df.loc[df['ice_thick'] > 0] - assert_allclose(np.abs(df['ice_flux'] - df['my_flux']), 0, atol=35e3) - assert df['ice_velocity'].max() > 25 - assert df['tributary_flux'].max() == 0 + # Also convert ours + df['ice_flux'] *= cfg.SEC_IN_YEAR + df['ice_velocity'] *= cfg.SEC_IN_YEAR + + assert_allclose(np.abs(df['ice_flux'] - df['my_flux']), 0, atol=35e3) + assert df['ice_velocity'].max() > 25 + + if isinstance(model, oggm.core.flowline.FluxBasedModel): + df['tributary_flux'] *= cfg.SEC_IN_YEAR + assert df['tributary_flux'].max() == 0 fls = dummy_width_bed_tributary() model = FluxBasedModel(fls, mb_model=mb, y0=0.) @@ -236,8 +256,9 @@ def test_min_slope(self): """ Check what is the min slope a flowline model can produce """ - models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel] - kwargs = [{'fixed_dt': 3 * SEC_IN_DAY}, {}, {}] + models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel, + SemiImplicitModel] + kwargs = [{'fixed_dt': 3 * SEC_IN_DAY}, {}, {}, {}] lens = [] surface_h = [] volume = [] @@ -270,49 +291,60 @@ def test_min_slope(self): surface_h.append(fls[-1].surface_h.copy()) np.testing.assert_allclose(lens[0][-1], lens[1][-1], atol=101) + np.testing.assert_allclose(lens[3][-1], lens[1][-1], atol=101) np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=2e-3) np.testing.assert_allclose(volume[1][-1], volume[2][-1], atol=5e-3) + np.testing.assert_allclose(volume[3][-1], volume[2][-1], atol=5e-3) assert utils.rmsd(volume[0], volume[2]) < 1e-2 assert utils.rmsd(volume[1], volume[2]) < 1e-2 + assert utils.rmsd(volume[3], volume[2]) < 1e-2 if do_plot: # pragma: no cover plt.figure() plt.plot(yrs, lens[0], 'r') plt.plot(yrs, lens[1], 'b') plt.plot(yrs, lens[2], 'g') + plt.plot(yrs, lens[3], 'm') plt.title('Compare Length') plt.xlabel('years') plt.ylabel('[m]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(yrs, volume[0], 'r') plt.plot(yrs, volume[1], 'b') plt.plot(yrs, volume[2], 'g') + plt.plot(yrs, volume[3], 'm') plt.title('Compare Volume') plt.xlabel('years') plt.ylabel('[km^3]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(yrs, min_slope[0], 'r') plt.plot(yrs, min_slope[1], 'b') plt.plot(yrs, min_slope[2], 'g') + plt.plot(yrs, min_slope[3], 'm') plt.title('Compare min slope') plt.xlabel('years') plt.ylabel('[degrees]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(fls[-1].bed_h, 'k') plt.plot(surface_h[0], 'r') plt.plot(surface_h[1], 'b') plt.plot(surface_h[2], 'g') + plt.plot(surface_h[3], 'm') plt.title('Compare Shape') plt.xlabel('[m]') plt.ylabel('Elevation [m]') - plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=3) + plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee', + 'Implicit'], loc=3) plt.show() @pytest.mark.slow @@ -322,7 +354,8 @@ def test_cliff(self): what the models do when the cliff height is changed """ - models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel] + models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel, + SemiImplicitModel] lens = [] surface_h = [] @@ -346,34 +379,40 @@ def test_cliff(self): volume.append(vol) surface_h.append(fls[-1].surface_h.copy()) - if False: # pragma: no cover + if do_plot: # pragma: no cover plt.figure() plt.plot(yrs, lens[0], 'r') plt.plot(yrs, lens[1], 'b') plt.plot(yrs, lens[2], 'g') + plt.plot(yrs, lens[3], 'm') plt.title('Compare Length') plt.xlabel('years') plt.ylabel('[m]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(yrs, volume[0], 'r') plt.plot(yrs, volume[1], 'b') plt.plot(yrs, volume[2], 'g') + plt.plot(yrs, volume[3], 'm') plt.title('Compare Volume') plt.xlabel('years') plt.ylabel('[km^3]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(fls[-1].bed_h, 'k') plt.plot(surface_h[0], 'r') plt.plot(surface_h[1], 'b') plt.plot(surface_h[2], 'g') + plt.plot(surface_h[3], 'm') plt.title('Compare Shape') plt.xlabel('[m]') plt.ylabel('Elevation [m]') - plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=3) + plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee', + 'Implicit'], loc=3) plt.show() # OK, so basically, Alex's tests below show that the other models @@ -384,6 +423,8 @@ def test_cliff(self): # Unit-testing perspective: # "verify" that indeed the models are wrong of more than 50% assert volume[1][-1] > volume[2][-1] * 1.5 + # SemiImplicit less wrong compared to FuxBased + assert volume[3][-1] > volume[2][-1] * 1.25 # Karthaus is even worse assert volume[0][-1] > volume[1][-1] @@ -403,7 +444,7 @@ def test_cliff(self): @pytest.mark.slow def test_equilibrium(self): - models = [KarthausModel, FluxBasedModel] + models = [KarthausModel, FluxBasedModel, SemiImplicitModel] vols = [] for model in models: @@ -433,8 +474,9 @@ def test_run_until(self): # Just check that exotic times are guaranteed to be met yrs = np.array([10.2, 10.2, 10.200001, 10.3, 99.999, 150.]) - models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel] - steps = [31 * SEC_IN_DAY, None, None] + models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel, + SemiImplicitModel] + steps = [31 * SEC_IN_DAY, None, None, None] # Annual update lens = [] @@ -462,13 +504,18 @@ def test_run_until(self): surface_h.append(fls[-1].surface_h.copy()) np.testing.assert_almost_equal(lens[0][-1], lens[1][-1]) + np.testing.assert_almost_equal(lens[3][-1], lens[1][-1]) np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2) np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2) + np.testing.assert_allclose(volume[3][-1], volume[2][-1], atol=1e-2) assert utils.rmsd(lens[0], lens[1]) < 50. + assert utils.rmsd(lens[0], lens[3]) < 50. assert utils.rmsd(volume[2], volume[1]) < 1e-3 + assert utils.rmsd(volume[2], volume[3]) < 1e-3 assert utils.rmsd(surface_h[0], surface_h[1]) < 5 assert utils.rmsd(surface_h[1], surface_h[2]) < 5 + assert utils.rmsd(surface_h[3], surface_h[2]) < 5 # Always update lens = [] @@ -493,19 +540,25 @@ def test_run_until(self): surface_h.append(fls[-1].surface_h.copy()) np.testing.assert_almost_equal(lens[0][-1], lens[1][-1]) + np.testing.assert_almost_equal(lens[0][-1], lens[3][-1]) np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2) np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2) + np.testing.assert_allclose(volume[0][-1], volume[3][-1], atol=1e-2) assert utils.rmsd(lens[0], lens[1]) < 50. + assert utils.rmsd(lens[0], lens[3]) < 50. assert utils.rmsd(volume[2], volume[1]) < 1e-3 + assert utils.rmsd(volume[2], volume[3]) < 1e-3 assert utils.rmsd(surface_h[0], surface_h[1]) < 5 assert utils.rmsd(surface_h[1], surface_h[2]) < 5 + assert utils.rmsd(surface_h[3], surface_h[2]) < 5 @pytest.mark.slow def test_adaptive_ts(self): - models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel] - steps = [31 * SEC_IN_DAY, None, None] + models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel, + SemiImplicitModel] + steps = [31 * SEC_IN_DAY, None, None, None] lens = [] surface_h = [] volume = [] @@ -528,19 +581,25 @@ def test_adaptive_ts(self): surface_h.append(fls[-1].surface_h.copy()) np.testing.assert_almost_equal(lens[0][-1], lens[1][-1]) + np.testing.assert_almost_equal(lens[0][-1], lens[3][-1]) np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2) np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2) + np.testing.assert_allclose(volume[0][-1], volume[3][-1], atol=1e-2) assert utils.rmsd(lens[0], lens[1]) < 50. + assert utils.rmsd(lens[0], lens[3]) < 50. assert utils.rmsd(volume[2], volume[1]) < 1e-3 + assert utils.rmsd(volume[2], volume[3]) < 1e-3 assert utils.rmsd(surface_h[0], surface_h[1]) < 5 assert utils.rmsd(surface_h[1], surface_h[2]) < 5 + assert utils.rmsd(surface_h[3], surface_h[2]) < 5 @pytest.mark.slow def test_bumpy_bed(self): - models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel] - steps = [15 * SEC_IN_DAY, None, None] + models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel, + SemiImplicitModel] + steps = [15 * SEC_IN_DAY, None, None, None] lens = [] surface_h = [] volume = [] @@ -567,46 +626,58 @@ def test_bumpy_bed(self): plt.plot(yrs, lens[0], 'r') plt.plot(yrs, lens[1], 'b') plt.plot(yrs, lens[2], 'g') + plt.plot(yrs, lens[3], 'm') plt.title('Compare Length') plt.xlabel('years') plt.ylabel('[m]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(yrs, volume[0], 'r') plt.plot(yrs, volume[1], 'b') plt.plot(yrs, volume[2], 'g') + plt.plot(yrs, volume[2], 'm') plt.title('Compare Volume') plt.xlabel('years') plt.ylabel('[km^3]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(fls[-1].bed_h, 'k') plt.plot(surface_h[0], 'r') plt.plot(surface_h[1], 'b') plt.plot(surface_h[2], 'g') + plt.plot(surface_h[2], 'm') plt.title('Compare Shape') plt.xlabel('[m]') plt.ylabel('Elevation [m]') - plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=3) + plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee', + 'Implicit'], loc=3) plt.show() np.testing.assert_almost_equal(lens[0][-1], lens[1][-1]) + np.testing.assert_almost_equal(lens[0][-1], lens[3][-1]) np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2) np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2) + np.testing.assert_allclose(volume[0][-1], volume[3][-1], atol=1e-2) assert utils.rmsd(lens[0], lens[1]) < 50. + assert utils.rmsd(lens[0], lens[3]) < 50. assert utils.rmsd(volume[0], volume[1]) < 1e-2 assert utils.rmsd(volume[0], volume[2]) < 1e-2 + assert utils.rmsd(volume[0], volume[3]) < 1e-2 assert utils.rmsd(surface_h[0], surface_h[1]) < 5 assert utils.rmsd(surface_h[0], surface_h[2]) < 5 + assert utils.rmsd(surface_h[0], surface_h[3]) < 5 @pytest.mark.slow def test_noisy_bed(self): - models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel] - steps = [15 * SEC_IN_DAY, None, None] + models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel, + SemiImplicitModel] + steps = [15 * SEC_IN_DAY, None, None, None] lens = [] surface_h = [] volume = [] @@ -634,40 +705,51 @@ def test_noisy_bed(self): plt.plot(yrs, lens[0], 'r') plt.plot(yrs, lens[1], 'b') plt.plot(yrs, lens[2], 'g') + plt.plot(yrs, lens[3], 'm') plt.title('Compare Length') plt.xlabel('years') plt.ylabel('[m]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(yrs, volume[0], 'r') plt.plot(yrs, volume[1], 'b') plt.plot(yrs, volume[2], 'g') + plt.plot(yrs, volume[2], 'm') plt.title('Compare Volume') plt.xlabel('years') plt.ylabel('[km^3]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(fls[-1].bed_h, 'k') plt.plot(surface_h[0], 'r') plt.plot(surface_h[1], 'b') plt.plot(surface_h[2], 'g') + plt.plot(surface_h[2], 'm') plt.title('Compare Shape') plt.xlabel('[m]') plt.ylabel('Elevation [m]') - plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=3) + plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee', + 'Implicit'], loc=3) plt.show() np.testing.assert_allclose(lens[0][-1], lens[1][-1], atol=101) + np.testing.assert_allclose(lens[0][-1], lens[3][-1], atol=101) np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2) np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2) + np.testing.assert_allclose(volume[0][-1], volume[3][-1], atol=1e-2) assert utils.rmsd(lens[0], lens[1]) < 100. + assert utils.rmsd(lens[0], lens[3]) < 100. assert utils.rmsd(volume[0], volume[1]) < 1e-1 assert utils.rmsd(volume[0], volume[2]) < 1e-1 + assert utils.rmsd(volume[0], volume[3]) < 1e-1 assert utils.rmsd(surface_h[0], surface_h[1]) < 10 assert utils.rmsd(surface_h[0], surface_h[2]) < 10 + assert utils.rmsd(surface_h[0], surface_h[3]) < 10 @pytest.mark.slow def test_varying_width(self): @@ -675,8 +757,9 @@ def test_varying_width(self): accumulation area twice as wide as the tongue.""" # set do_plot = True to see the plots - models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel] - steps = [15 * SEC_IN_DAY, None, None] + models = [KarthausModel, FluxBasedModel, MUSCLSuperBeeModel, + SemiImplicitModel] + steps = [15 * SEC_IN_DAY, None, None, None] lens = [] surface_h = [] volume = [] @@ -703,39 +786,52 @@ def test_varying_width(self): plt.plot(yrs, lens[0], 'r') plt.plot(yrs, lens[1], 'b') plt.plot(yrs, lens[2], 'g') + plt.plot(yrs, lens[3], 'm') plt.title('Compare Length') plt.xlabel('years') plt.ylabel('[m]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(yrs, volume[0], 'r') plt.plot(yrs, volume[1], 'b') plt.plot(yrs, volume[2], 'g') + plt.plot(yrs, volume[3], 'm') plt.title('Compare Volume') plt.xlabel('years') plt.ylabel('[km^3]') - plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=2) + plt.legend(['Karthaus', 'Flux', 'MUSCL-SuperBee', 'Implicit'], + loc=2) plt.figure() plt.plot(fls[-1].bed_h, 'k') plt.plot(surface_h[0], 'r') plt.plot(surface_h[1], 'b') plt.plot(surface_h[2], 'g') + plt.plot(surface_h[3], 'm') plt.title('Compare Shape') plt.xlabel('[m]') plt.ylabel('Elevation [m]') - plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee'], loc=3) + plt.legend(['Bed', 'Karthaus', 'Flux', 'MUSCL-SuperBee', + 'Implicit'], loc=3) plt.show() np.testing.assert_almost_equal(lens[0][-1], lens[1][-1]) + np.testing.assert_almost_equal(lens[0][-1], lens[3][-1]) np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=2e-2) + np.testing.assert_allclose(volume[1][-1], volume[3][-1], atol=1e-2) np.testing.assert_allclose(utils.rmsd(lens[0], lens[1]), 0., atol=70) + np.testing.assert_allclose(utils.rmsd(lens[3], lens[1]), 0., atol=50) np.testing.assert_allclose(utils.rmsd(volume[0], volume[1]), 0., atol=1e-2) + np.testing.assert_allclose(utils.rmsd(volume[3], volume[1]), 0., + atol=4e-3) np.testing.assert_allclose(utils.rmsd(surface_h[0], surface_h[1]), 0., atol=5) + np.testing.assert_allclose(utils.rmsd(surface_h[3], surface_h[1]), 0., + atol=3) @pytest.mark.slow def test_tributary(self): @@ -860,8 +956,9 @@ def test_trapezoidal_bed(self): akm = (tb._w0_m + tb._lambdas * h) * len(sec) * 100 np.testing.assert_almost_equal(tb.area_m2, akm) - models = [KarthausModel, FluxBasedModel] - flss = [dummy_constant_bed(), dummy_trapezoidal_bed()] + models = [KarthausModel, FluxBasedModel, SemiImplicitModel] + flss = [dummy_constant_bed(), dummy_trapezoidal_bed(), + dummy_trapezoidal_bed()] lens = [] surface_h = [] @@ -888,23 +985,41 @@ def test_trapezoidal_bed(self): surface_h.append(fls[-1].surface_h.copy()) np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=1e-2) + np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=1e-2) + + np.testing.assert_allclose(lens[2][-1], lens[1][-1]) + np.testing.assert_allclose(volume[2][-1], volume[1][-1], atol=2e-5) + + np.testing.assert_allclose(utils.rmsd(lens[2], lens[1]), 0., atol=6) + np.testing.assert_allclose(utils.rmsd(volume[2], volume[1]), 0., + atol=5e-5) + np.testing.assert_allclose(utils.rmsd(surface_h[2], surface_h[1]), 0., + atol=2e-2) if do_plot: # pragma: no cover plt.plot(lens[0], 'r') plt.plot(lens[1], 'b') + plt.plot(lens[2], 'm') + plt.title('Length') plt.show() plt.plot(volume[0], 'r') plt.plot(volume[1], 'b') + plt.plot(volume[2], 'm') + plt.title('Volume') plt.show() plt.plot(fls[-1].bed_h, 'k') plt.plot(surface_h[0], 'r') plt.plot(surface_h[1], 'b') + plt.plot(surface_h[2], 'm') + plt.title('Surface_h') plt.show() plt.plot(widths[0], 'r') plt.plot(widths[1], 'b') + plt.plot(widths[2], 'm') + plt.title('Widths') plt.show() @pytest.mark.slow @@ -1011,6 +1126,80 @@ def test_mixed_bed(self): plt.legend() plt.show() + @pytest.mark.slow + def test_mixed_trap_rect_bed(self): + + models = [KarthausModel, FluxBasedModel, SemiImplicitModel] + flss = [dummy_constant_bed(), dummy_mixed_trap_rect_bed(), + dummy_mixed_trap_rect_bed()] + + lens = [] + surface_h = [] + volume = [] + widths = [] + yrs = np.arange(1, 700, 2) + for model, fls in zip(models, flss): + mb = LinearMassBalance(2800.) + + model = model(fls, mb_model=mb, fs=self.fs_old, + glen_a=self.aglen_old, + fixed_dt=14 * SEC_IN_DAY) + + length = yrs * 0. + vol = yrs * 0. + for i, y in enumerate(yrs): + model.run_until(y) + assert model.yr == y + length[i] = fls[-1].length_m + vol[i] = fls[-1].volume_km3 + lens.append(length) + volume.append(vol) + widths.append(fls[-1].widths_m.copy()) + surface_h.append(fls[-1].surface_h.copy()) + + np.testing.assert_allclose(volume[0][-1], volume[1][-1], atol=2e-1) + np.testing.assert_allclose(volume[0][-1], volume[2][-1], atol=2e-1) + + # maximum allow difference in length is one grid point + np.testing.assert_allclose(lens[2][-1], lens[1][-1], atol=1e2) + + np.testing.assert_allclose(volume[2][-1], volume[1][-1], atol=2e-5) + np.testing.assert_allclose(utils.rmsd(lens[2], lens[1]), 0., atol=23) + np.testing.assert_allclose(utils.rmsd(volume[2], volume[1]), 0., + atol=3e-5) + np.testing.assert_allclose(utils.rmsd(surface_h[2], surface_h[1]), 0., + atol=6e-2) + + if do_plot: # pragma: no cover + plt.plot(lens[0], 'r', label='normal') + plt.plot(lens[1], 'b', label='mixed flux') + plt.plot(lens[2], 'm', label='mixed impl') + plt.title('Length') + plt.legend() + plt.show() + + plt.plot(volume[0], 'r', label='normal') + plt.plot(volume[1], 'b', label='mixed flux') + plt.plot(volume[2], 'm', label='mixed impl') + plt.title('Volume') + plt.legend() + plt.show() + + plt.plot(fls[-1].bed_h, 'k') + plt.plot(surface_h[0], 'r', label='normal') + plt.plot(surface_h[1], 'b', label='mixed flux') + plt.plot(surface_h[2], 'm', label='mixed impl') + plt.title('Surface_h') + plt.legend() + plt.show() + + plt.plot(widths[0], 'r', label='normal') + plt.plot(widths[1], 'b', label='mixed flux') + plt.plot(widths[2], 'm', label='mixed impl') + plt.title('Widths') + plt.legend() + plt.show() + @pytest.mark.slow def test_raise_on_boundary(self): @@ -1316,18 +1505,22 @@ def setUp(self): def test_find_flux_from_thickness(self): mb = LinearMassBalance(2600.) - model = FluxBasedModel(dummy_constant_bed(), mb_model=mb) - model.run_until(700) - - # Pick a flux and slope somewhere in the glacier - for i in [1, 10, 20, 50]: - flux = model.flux_stag[0][i] - slope = model.slope_stag[0][i] - thick = model.thick_stag[0][i] - width = model.fls[0].widths_m[i] - - out = find_sia_flux_from_thickness(slope, width, thick) - assert_allclose(out, flux, atol=1e-7) + model_flux = FluxBasedModel(dummy_constant_bed(), mb_model=mb) + model_flux.run_until(700) + + model_impl = SemiImplicitModel(dummy_constant_bed(), mb_model=mb) + model_impl.run_until(700) + + for model in [model_flux, model_impl]: + # Pick a flux and slope somewhere in the glacier + for i in [1, 10, 20, 50]: + flux = model.flux_stag[0][i] + slope = model.slope_stag[0][i] + thick = model.thick_stag[0][i] + width = model.fls[0].widths_m[i] + + out = find_sia_flux_from_thickness(slope, width, thick) + assert_allclose(out, flux, atol=1e-7) def test_simple_flux_gate(self):