diff --git a/test/tracerEq/test_consistency_2d.py b/test/tracerEq/test_consistency_2d.py new file mode 100644 index 000000000..28595e0a2 --- /dev/null +++ b/test/tracerEq/test_consistency_2d.py @@ -0,0 +1,128 @@ +# Tracer box in 2D +# ================ +# +# Solves a standing wave in a rectangular basin using wave equation. +# +# This version uses a constant tracer to check local/global conservation of tracers. +# +# Initial condition for elevation corresponds to a standing wave. +# Time step and export interval are chosen based on theorethical +# oscillation frequency. Initial condition repeats every 20 exports. +# +# +# Thanasis Angeloudis based on Tuomas Karna test_consistency.py 15-01-2018 + +from thetis import * + +def run_tracer_consistency(constant_c = True, **model_options): + + t_cycle = 2000.0 # standing wave period + depth = 50.0 # average depth + lx = np.sqrt(9.81*depth)*t_cycle # wave length + ly = 3000.0 + nx = 18 + ny = 2 + mesh2d = RectangleMesh(nx, ny, lx, ly) + tracer_value = 4.5 + elev_amp = 2.0 + # estimate of max advective velocity used to estimate time step + u_mag = Constant(1.0) + + outputdir = 'outputs' + + # bathymetry + p1_2d = FunctionSpace(mesh2d, 'CG', 1) + bathymetry_2d = Function(p1_2d, name='Bathymetry') + x_2d, y_2d = SpatialCoordinate(mesh2d) + # non-trivial bathymetry, to properly test 2d tracer conservation + bathymetry_2d.interpolate(depth + depth/10.*sin(x_2d/lx*pi)) + + + # set time step, export interval and run duration + n_steps = 8 + t_export = round(float(t_cycle/n_steps)) + # for testing tracer conservation, we don't want to come back to the initial condition + t_end = 2.5*t_cycle + + # create solver + solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d) + options = solver_obj.options + options.use_nonlinear_equations = True + options.solve_tracer = True + options.use_limiter_for_tracers = not constant_c + options.simulation_export_time = t_export + options.simulation_end_time = t_end + options.horizontal_velocity_scale = Constant(u_mag) + options.check_volume_conservation_2d = True + options.check_tracer_conservation = True + options.check_tracer_overshoot = True + options.timestepper_type = 'CrankNicolson' + options.output_directory = outputdir + options.fields_to_export = ['uv_2d', 'elev_2d', 'tracer_2d'] + options.update(model_options) + + if not options.no_exports: + print_output('Exporting to {:}'.format(options.output_directory)) + + solver_obj.create_function_spaces() + elev_init = Function(solver_obj.function_spaces.H_2d) + elev_init.project(-elev_amp*cos(2*pi*x_2d/lx)) + + tracer_init2d = None + if options.solve_tracer and constant_c: + tracer_init2d = Function(solver_obj.function_spaces.Q_2d, name='initial tracer') + tracer_init2d.assign(tracer_value) + if options.solve_tracer and not constant_c: + tracer_init2d = Function(solver_obj.function_spaces.Q_2d, name='initial tracer') + tracer_l = 0 + tracer_r = 30.0 + tracer_init2d.interpolate(tracer_l + (tracer_r - tracer_l)*0.5*(1.0 + sign(x_2d - lx/4))) + + solver_obj.assign_initial_conditions(elev=elev_init, tracer=tracer_init2d ) + solver_obj.iterate() + + # TODO do these checks every export ... + vol2d, vol2d_rerr = solver_obj.callbacks['export']['volume2d']() + assert vol2d_rerr < 1e-10, '2D volume is not conserved' + if options.solve_tracer: + tracer_int, tracer_int_rerr = solver_obj.callbacks['export']['tracer_2d mass']() + assert abs(tracer_int_rerr) < 1e-4, 'tracer is not conserved' + smin, smax, undershoot, overshoot = solver_obj.callbacks['export']['tracer_2d overshoot']() + max_abs_overshoot = max(abs(undershoot), abs(overshoot)) + overshoot_tol = 1e-12 + msg = 'Tracer overshoots are too large: {:}'.format(max_abs_overshoot) + assert max_abs_overshoot < overshoot_tol, msg + + + +def test_const_tracer(): + """ + Test CrankNicolson timeintegrator without slope limiters + Constant tracer, should remain constant + """ + run_tracer_consistency(constant_c= True, + use_nonlinear_equations=True, + solve_tracer=True, + use_limiter_for_tracers=False, + no_exports=True) + + + +def test_nonconst_tracer(): + """ + Test CrankNicolson timeintegrator with slope limiters + Non-trivial tracer, should see no overshoots and be conserved + """ + run_tracer_consistency(constant_c= False, + use_nonlinear_equations=True, + solve_tracer=True, + use_limiter_for_tracers=True, + no_exports=True) + + +if __name__ == '__main__': + run_tracer_consistency(constant_c= False, + use_nonlinear_equations=True, + solve_tracer=True, + use_limiter_for_tracers=False, + no_exports=False) diff --git a/test/tracerEq/test_h-advection_mes.py b/test/tracerEq/test_h-advection_mes.py index bbe7e954e..5ca95ebb8 100644 --- a/test/tracerEq/test_h-advection_mes.py +++ b/test/tracerEq/test_h-advection_mes.py @@ -120,7 +120,7 @@ def export_func(): next_export_t += solverobj.options.simulation_export_time iexport += 1 - # project analytical solultion on high order mesh + # project analytical solution on high order mesh t_const.assign(t) salt_ana_ho.project(ana_salt_expr) # compute L2 norm diff --git a/test/tracerEq/test_h-advection_mes_2d.py b/test/tracerEq/test_h-advection_mes_2d.py new file mode 100644 index 000000000..8de74e3a8 --- /dev/null +++ b/test/tracerEq/test_h-advection_mes_2d.py @@ -0,0 +1,198 @@ +""" +Testing 2D horizontal advection of tracers + +Tuomas Karna, edited by Athanasios Angeloudis +""" +from thetis import * +import numpy +from scipy import stats +import pytest + + +def run(refinement, **model_options): + print_output('--- running refinement {:}'.format(refinement)) + + # domain dimensions - channel in x-direction + lx = 15.0e3 + ly = 6.0e3/refinement + area = lx*ly + depth = 40.0 + u = 1.0 + + # mesh + nx = 6*refinement + 1 + ny = 1 # constant -- channel + mesh2d = RectangleMesh(nx, ny, lx, ly) + + # simulation run time + t_end = 3000.0 + # initial time + t_init = 0.0 + t_export = (t_end - t_init)/8.0 + + # outputs + outputdir = 'outputs' + + # bathymetry + P1_2d = FunctionSpace(mesh2d, 'CG', 1) + bathymetry_2d = Function(P1_2d, name='Bathymetry') + bathymetry_2d.assign(depth) + + solverobj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d) + options = solverobj.options + options.use_nonlinear_equations = False + options.use_lax_friedrichs_velocity = True + options.lax_friedrichs_velocity_scaling_factor = Constant(1.0) + options.use_lax_friedrichs_tracer = False + options.horizontal_velocity_scale = Constant(abs(u)) + options.no_exports = True + options.output_directory = outputdir + options.simulation_end_time = t_end + options.simulation_export_time = t_export + options.solve_tracer = True + options.use_limiter_for_tracers = True + options.fields_to_export = ['tracer_2d'] + options.update(model_options) + uv_expr = as_vector((u, 0)) + bnd_salt_2d = {'value': Constant(0.0), 'uv': uv_expr} + bnd_uv_2d = {'uv': uv_expr} + solverobj.bnd_functions['tracer'] = { + 1: bnd_salt_2d, + 2: bnd_salt_2d, + } + solverobj.bnd_functions['momentum'] = { + 1: bnd_uv_2d, + 2: bnd_uv_2d, + } + + solverobj.create_equations() + + t = t_init # simulation time + + x0 = 0.3*lx + sigma = 1600. + xy = SpatialCoordinate(solverobj.mesh2d) + t_const = Constant(t) + ana_tracer_expr = exp(-(xy[0] - x0 - u*t_const)**2/sigma**2) + tracer_ana = Function(solverobj.function_spaces.H_2d, name='tracer analytical') + uv_init = Function(solverobj.function_spaces.U_2d, name='initial uv') + uv_init.project(uv_expr) + solverobj.assign_initial_conditions(uv=uv_init, tracer=ana_tracer_expr) + + # export analytical solution + if not options.no_exports: + out_tracer_ana = File(os.path.join(options.output_directory, 'tracer_ana.pvd')) + + def export_func(): + if not options.no_exports: + solverobj.export() + # update analytical solution to correct time + t_const.assign(t) + out_tracer_ana.write(tracer_ana.project(ana_tracer_expr)) + + # export initial conditions + export_func() + + # custom time loop that solves tracer equation only + ti = solverobj.timestepper.timesteppers.tracer + + i = 0 + iexport = 1 + next_export_t = t + solverobj.options.simulation_export_time + while t < t_end - 1e-8: + ti.advance(t) + t += solverobj.dt + i += 1 + if t >= next_export_t - 1e-8: + print_output('{:3d} i={:5d} t={:8.2f} s salt={:8.2f}'.format(iexport, i, t, norm(solverobj.fields.tracer_2d))) + export_func() + next_export_t += solverobj.options.simulation_export_time + iexport += 1 + + # project analytical solution on high order mesh + t_const.assign(t) + + # compute L2 norm + l2_err = errornorm(ana_tracer_expr, solverobj.fields.tracer_2d)/numpy.sqrt(area) + print_output('L2 error {:.12f}'.format(l2_err)) + + return l2_err + + +def run_convergence(ref_list, saveplot=False, **options): + """Runs test for a list of refinements and computes error convergence rate""" + + polynomial_degree = options.get('polynomial_degree', 1) + l2_err = [] + for r in ref_list: + l2_err.append(run(r, **options)) + x_log = numpy.log10(numpy.array(ref_list, dtype=float)**-1) + y_log = numpy.log10(numpy.array(l2_err)) + setup_name = 'h-diffusion' + + def check_convergence(x_log, y_log, expected_slope, field_str, saveplot): + slope_rtol = 0.20 + slope, intercept, r_value, p_value, std_err = stats.linregress(x_log, y_log) + if saveplot: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(figsize=(5, 5)) + # plot points + ax.plot(x_log, y_log, 'k.') + x_min = x_log.min() + x_max = x_log.max() + offset = 0.05*(x_max - x_min) + npoints = 50 + xx = numpy.linspace(x_min - offset, x_max + offset, npoints) + yy = intercept + slope*xx + # plot line + ax.plot(xx, yy, linestyle='--', linewidth=0.5, color='k') + ax.text(xx[2*int(npoints/3)], yy[2*int(npoints/3)], '{:4.2f}'.format(slope), + verticalalignment='top', + horizontalalignment='left') + ax.set_xlabel('log10(dx)') + ax.set_ylabel('log10(L2 error)') + ax.set_title(' '.join([setup_name, field_str, 'degree={:}'.format(polynomial_degree)])) + ref_str = 'ref-' + '-'.join([str(r) for r in ref_list]) + degree_str = 'o{:}'.format(polynomial_degree) + imgfile = '_'.join(['convergence', setup_name, field_str, ref_str, degree_str]) + imgfile += '.png' + imgdir = create_directory('plots') + imgfile = os.path.join(imgdir, imgfile) + print_output('saving figure {:}'.format(imgfile)) + plt.savefig(imgfile, dpi=200, bbox_inches='tight') + if expected_slope is not None: + err_msg = '{:}: Wrong convergence rate {:.4f}, expected {:.4f}'.format(setup_name, slope, expected_slope) + assert slope > expected_slope*(1 - slope_rtol), err_msg + print_output('{:}: convergence rate {:.4f} PASSED'.format(setup_name, slope)) + else: + print_output('{:}: {:} convergence rate {:.4f}'.format(setup_name, field_str, slope)) + return slope + + check_convergence(x_log, y_log, polynomial_degree+1, 'tracer', saveplot) + +# --------------------------- +# standard tests for pytest +# --------------------------- + + +@pytest.fixture(params=[1]) +def polynomial_degree(request): + return request.param + + +@pytest.mark.parametrize(('stepper'), + [('CrankNicolson')]) + +def test_horizontal_advection(polynomial_degree, stepper, ): + run_convergence([1, 2, 3], polynomial_degree=polynomial_degree, + timestepper_type=stepper) + +# --------------------------- +# run individual setup for debugging +# --------------------------- + + +if __name__ == '__main__': + run_convergence([1, 2, 3, 4], polynomial_degree=1, + timestepper_type='CrankNicolson', + no_exports=False, saveplot=True) diff --git a/test/tracerEq/test_h-diffusion_mes_2d.py b/test/tracerEq/test_h-diffusion_mes_2d.py new file mode 100644 index 000000000..66b9ee4f8 --- /dev/null +++ b/test/tracerEq/test_h-diffusion_mes_2d.py @@ -0,0 +1,176 @@ +""" +Testing 2D horizontal diffusion of tracers against analytical solution. + +""" +from thetis import * +import numpy +from scipy import stats +import pytest + + +def run(refinement, **model_options): + print_output('--- running refinement {:}'.format(refinement)) + + # domain dimensions - channel in x-direction + lx = 20.0e3 + ly = 5.0e3 / refinement + area = lx*ly + depth = 30.0 + horizontal_diffusivity = 1.0e3 + + # mesh + nx = 8 * refinement + 1 + ny = 1 # constant -- channel + mesh2d = RectangleMesh(nx, ny, lx, ly) + + # simulation run time + t_end = 3000.0 + # initial time + t_init = 1000.0 # NOTE start from t > 0 for smoother init cond + t_export = (t_end - t_init)/8.0 + + # outputs + outputdir = 'outputs' + + # bathymetry + p1_2d = FunctionSpace(mesh2d, 'CG', 1) + bathymetry_2d = Function(p1_2d, name='Bathymetry') + bathymetry_2d.assign(depth) + + solverobj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d) + options = solverobj.options + options.use_nonlinear_equations = False + options.horizontal_velocity_scale = Constant(1.0) + options.no_exports = True + options.output_directory = outputdir + options.simulation_end_time = t_end + options.simulation_export_time = t_export + options.solve_tracer = True + options.use_limiter_for_tracers = True + options.fields_to_export = ['tracer_2d'] + options.horizontal_diffusivity = Constant(horizontal_diffusivity) + options.horizontal_viscosity_scale = Constant(horizontal_diffusivity) + options.update(model_options) + + solverobj.create_equations() + + t = t_init # simulation time + + t_const = Constant(t) + u_max = 1.0 + u_min = -1.0 + x0 = lx/2.0 + x, y = SpatialCoordinate(solverobj.mesh2d) + tracer_expr = 0.5*(u_max + u_min) - 0.5*(u_max - u_min)*erf((x - x0)/sqrt(4*horizontal_diffusivity*t_const)) + tracer_ana = Function(solverobj.function_spaces.H_2d, name='tracer analytical') + elev_init = Function(solverobj.function_spaces.H_2d, name='elev init') + solverobj.assign_initial_conditions(elev=elev_init, tracer=tracer_expr) + + # export analytical solution + if not options.no_exports: + out_tracer_ana = File(os.path.join(options.output_directory, 'tracer_ana.pvd')) + + def export_func(): + if not options.no_exports: + solverobj.export() + # update analytical solution to correct time + t_const.assign(t) + out_tracer_ana.write(tracer_ana.project(tracer_expr)) + + # export initial conditions + export_func() + + # custom time loop that solves tracer equation only + ti = solverobj.timestepper.timesteppers.tracer + + i = 0 + iexport = 1 + next_export_t = t + solverobj.options.simulation_export_time + while t < t_end - 1e-8: + ti.advance(t) + t += solverobj.dt + i += 1 + if t >= next_export_t - 1e-8: + print_output('{:3d} i={:5d} t={:8.2f} s tracer={:8.2f}'.format(iexport, i, t, norm(solverobj.fields.tracer_2d))) + export_func() + next_export_t += solverobj.options.simulation_export_time + iexport += 1 + + # project analytical solution on high order mesh + t_const.assign(t) + # compute L2 norm + l2_err = errornorm(tracer_expr, solverobj.fields.tracer_2d)/numpy.sqrt(area) + print_output('L2 error {:.12f}'.format(l2_err)) + + return l2_err + + +def run_convergence(ref_list, saveplot=False, **options): + """Runs test for a list of refinements and computes error convergence rate""" + polynomial_degree = options.get('polynomial_degree', 1) + l2_err = [] + for r in ref_list: + l2_err.append(run(r, **options)) + x_log = numpy.log10(numpy.array(ref_list, dtype=float)**-1) + y_log = numpy.log10(numpy.array(l2_err)) + setup_name = 'h-diffusion' + + def check_convergence(x_log, y_log, expected_slope, field_str, saveplot): + slope_rtol = 0.20 + slope, intercept, r_value, p_value, std_err = stats.linregress(x_log, y_log) + if saveplot: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(figsize=(5, 5)) + # plot points + ax.plot(x_log, y_log, 'k.') + x_min = x_log.min() + x_max = x_log.max() + offset = 0.05*(x_max - x_min) + npoints = 50 + xx = numpy.linspace(x_min - offset, x_max + offset, npoints) + yy = intercept + slope*xx + # plot line + ax.plot(xx, yy, linestyle='--', linewidth=0.5, color='k') + ax.text(xx[2*int(npoints/3)], yy[2*int(npoints/3)], '{:4.2f}'.format(slope), + verticalalignment='top', + horizontalalignment='left') + ax.set_xlabel('log10(dx)') + ax.set_ylabel('log10(L2 error)') + ax.set_title(' '.join([setup_name, field_str, 'degree={:}'.format(polynomial_degree)])) + ref_str = 'ref-' + '-'.join([str(r) for r in ref_list]) + degree_str = 'o{:}'.format(polynomial_degree) + imgfile = '_'.join(['convergence', setup_name, field_str, ref_str, degree_str]) + imgfile += '.png' + imgdir = create_directory('plots') + imgfile = os.path.join(imgdir, imgfile) + print_output('saving figure {:}'.format(imgfile)) + plt.savefig(imgfile, dpi=200, bbox_inches='tight') + if expected_slope is not None: + err_msg = '{:}: Wrong convergence rate {:.4f}, expected {:.4f}'.format(setup_name, slope, expected_slope) + assert slope > expected_slope*(1 - slope_rtol), err_msg + print_output('{:}: convergence rate {:.4f} PASSED'.format(setup_name, slope)) + else: + print_output('{:}: {:} convergence rate {:.4f}'.format(setup_name, field_str, slope)) + return slope + + check_convergence(x_log, y_log, polynomial_degree+1, 'tracer', saveplot) + +# --------------------------- +# standard tests for pytest +# --------------------------- + +@pytest.mark.parametrize(('stepper'), + [('CrankNicolson')]) +def test_horizontal_diffusion(stepper): + run_convergence([1, 2, 3], polynomial_degree=1, + timestepper_type=stepper, + ) +# --------------------------- +# run individual setup for debugging +# --------------------------- + + +if __name__ == '__main__': + run_convergence([1, 2, 4, 6, 8], polynomial_degree=1, + timestepper_type='CrankNicolson', + no_exports=False, saveplot=True) diff --git a/test/tracerEq/test_steady_adv-diff_mms_2d.py b/test/tracerEq/test_steady_adv-diff_mms_2d.py new file mode 100644 index 000000000..0b2774170 --- /dev/null +++ b/test/tracerEq/test_steady_adv-diff_mms_2d.py @@ -0,0 +1,215 @@ +""" +Testing 2D tracer advection-diffusion equation with method of manufactured solution (MMS). + +Tuomas Karna 2015-11-28 +""" +from thetis import * +import numpy +from scipy import stats +import pytest + + +class Setup1: + """ + Constant bathymetry and u velocty, zero diffusivity, non-trivial tracer + """ + def bath(self, x, y, lx, ly): + return Constant(40.0) + + def elev(self, x, y, lx, ly): + return Constant(0.0) + + def uv(self, x, y, lx, ly): + return as_vector( + [ + Constant(1.0), + Constant(0.0), + ]) + + def kappa(self, x, y, lx, ly): + return Constant(0.0) + + def tracer(self, x, y, lx, ly): + return sin(0.2*pi*(3.0*x + 1.0*y)/lx) + + def residual(self, x, y, lx, ly): + return 0.6*pi*cos(0.2*pi*(3.0*x + 1.0*y)/lx)/lx + + +class Setup2: + """ + Constant bathymetry, zero velocity, constant kappa, x-varying T + """ + def bath(self, x, y, lx, ly): + return Constant(40.0) + + def elev(self, x, y, lx, ly): + return Constant(0.0) + + def uv(self, x, y, lx, ly): + return as_vector( + [ + Constant(1.0), + Constant(0.0), + ]) + + def kappa(self, x, y, lx, ly): + return Constant(50.0) + + def tracer(self, x, y, lx, ly): + return sin(3*pi*x/lx) + + def residual(self, x, y, lx, ly): + return 3.0*pi*cos(3*pi*x/lx)/lx - 450.0*pi**2*sin(3*pi*x/lx)/lx**2 + + +def run(setup, refinement, do_export=True, **options): + """Run single test and return L2 error""" + print_output('--- running {:} refinement {:}'.format(setup.__name__, refinement)) + # domain dimensions + lx = 15e3 + ly = 10e3 + area = lx*ly + t_end = 200.0 + + setup_obj = setup() + + # mesh + nx = 4*refinement + ny = 4*refinement + mesh2d = RectangleMesh(nx, ny, lx, ly) + + # outputs + outputdir = 'outputs' + if do_export: + out_2 = File(os.path.join(outputdir, 'Tracer.pvd')) + + # bathymetry + x_2d, y_2d = SpatialCoordinate(mesh2d) + p1_2d = FunctionSpace(mesh2d, 'CG', 1) + bathymetry_2d = Function(p1_2d, name='Bathymetry') + bathymetry_2d.project(setup_obj.bath(x_2d, y_2d, lx, ly)) + + solver_obj = solver2d.FlowSolver2d(mesh2d, bathymetry_2d, ) + solver_obj.options.element_family = 'dg-dg' + solver_obj.options.horizontal_velocity_scale = Constant(1.0) + solver_obj.options.no_exports = not do_export + solver_obj.options.output_directory = outputdir + solver_obj.options.simulation_end_time = t_end + solver_obj.options.fields_to_export = ['tracer_2d', 'uv_2d',] + solver_obj.options.horizontal_viscosity_scale = Constant(50.0) + solver_obj.options.update(options) + solver_obj.options.solve_tracer = True + solver_obj.options.timestepper_options.implicitness_theta = 1.0 + solver_obj.create_function_spaces() + + # functions for source terms + x, y= SpatialCoordinate(solver_obj.mesh2d) + solver_obj.options.tracer_source_2d = setup_obj.residual(x, y, lx, ly) + + # diffusivuty + solver_obj.options.horizontal_diffusivity = setup_obj.kappa(x, y, lx, ly) + + # analytical solution + trac_ana = setup_obj.tracer(x, y, lx, ly) + + bnd_tracer = {'value': trac_ana} + solver_obj.bnd_functions['tracer'] = {1: bnd_tracer, 2: bnd_tracer, + 3: bnd_tracer, 4: bnd_tracer} + + solver_obj.create_equations() + solver_obj.assign_initial_conditions(elev = setup_obj.elev(x, y, lx, ly), + uv = setup_obj.uv(x, y, lx, ly), + tracer = setup_obj.tracer(x, y, lx, ly)) + + # solve tracer advection-diffusion equation with residual source term + ti = solver_obj.timestepper + ti.timesteppers.tracer.initialize(ti.fields.tracer_2d) + + t = 0 + while t < t_end : + ti.timesteppers.tracer.advance(t) + if ti.options.use_limiter_for_tracers: + ti.solver.tracer_limiter.apply(ti.fields.tracer_2d) + t += solver_obj.dt + if do_export: + out_2.write(ti.fields.tracer_2d) + + l2_err = errornorm(trac_ana, solver_obj.fields.tracer_2d)/numpy.sqrt(area) + print_output('L2 error {:.12f}'.format(l2_err)) + + return l2_err + + +def run_convergence(setup, ref_list, do_export=False, save_plot=False, **options): + """Runs test for a list of refinements and computes error convergence rate""" + l2_err = [] + for r in ref_list: + l2_err.append(run(setup, r, do_export=do_export, **options)) + x_log = numpy.log10(numpy.array(ref_list, dtype=float)**-1) + y_log = numpy.log10(numpy.array(l2_err)) + + def check_convergence(x_log, y_log, expected_slope, field_str, save_plot): + slope_rtol = 0.2 + slope, intercept, r_value, p_value, std_err = stats.linregress(x_log, y_log) + setup_name = setup.__name__ + if save_plot: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(figsize=(5, 5)) + # plot points + ax.plot(x_log, y_log, 'k.') + x_min = x_log.min() + x_max = x_log.max() + offset = 0.05*(x_max - x_min) + npoints = 50 + xx = numpy.linspace(x_min - offset, x_max + offset, npoints) + yy = intercept + slope*xx + # plot line + ax.plot(xx, yy, linestyle='--', linewidth=0.5, color='k') + ax.text(xx[2*int(npoints/3)], yy[2*int(npoints/3)], '{:4.2f}'.format(slope), + verticalalignment='top', + horizontalalignment='left') + ax.set_xlabel('log10(dx)') + ax.set_ylabel('log10(L2 error)') + ax.set_title('tracer adv-diff MMS DG ') + ref_str = 'ref-' + '-'.join([str(r) for r in ref_list]) + + imgfile = '_'.join(['convergence', setup_name, field_str, ref_str,]) + imgfile += '.png' + img_dir = create_directory('plots') + imgfile = os.path.join(img_dir, imgfile) + print_output('saving figure {:}'.format(imgfile)) + plt.savefig(imgfile, dpi=200, bbox_inches='tight') + if expected_slope is not None: + err_msg = '{:}: Wrong convergence rate {:.4f}, expected {:.4f}'.format(setup_name, slope, expected_slope) + assert abs(slope - expected_slope)/expected_slope < slope_rtol, err_msg + print_output('{:}: convergence rate {:.4f} PASSED'.format(setup_name, slope)) + else: + print_output('{:}: {:} convergence rate {:.4f}'.format(setup_name, field_str, slope)) + return slope + + check_convergence(x_log, y_log, 2, 'tracer', save_plot) + + +# --------------------------- +# standard tests for pytest +# --------------------------- + +@pytest.fixture(params=['CrankNicolson']) +def timestepper_type(request): + return request.param + + +@pytest.fixture(params=[Setup1, + Setup2,], + ids=['setup1', 'setup2',]) +def setup(request): + return request.param + + +def test_convergence(setup, timestepper_type): + run_convergence(setup, [1, 2, 3], save_plot=False, timestepper_type=timestepper_type) + + +if __name__ == '__main__': + run_convergence(Setup1, [1, 2, 3], save_plot=True, timestepper_type='CrankNicolson') diff --git a/thetis/__init__.py b/thetis/__init__.py index 5e3e771c0..b528fcb11 100644 --- a/thetis/__init__.py +++ b/thetis/__init__.py @@ -11,6 +11,7 @@ import thetis.timezone as timezone # NOQA from thetis._version import get_versions from thetis.assembledschur import AssembledSchurPC # NOQA + __version__ = get_versions()['version'] del get_versions diff --git a/thetis/callback.py b/thetis/callback.py index 1602c2c6c..8cdad4089 100644 --- a/thetis/callback.py +++ b/thetis/callback.py @@ -301,12 +301,32 @@ def __init__(self, solver_obj, **kwargs): :arg solver_obj: Thetis solver object :arg **kwargs: any additional keyword arguments, see DiagnosticCallback """ + def vol2d(): return comp_volume_2d(self.solver_obj.fields.elev_2d, self.solver_obj.fields.bathymetry_2d) super(VolumeConservation2DCallback, self).__init__(vol2d, solver_obj, **kwargs) +class TracerMassConservation2DCallback(ScalarConservationCallback): + """Checks conservation of depth-averaged tracer mass (integral of 2D tracer multiplied by total depth)""" + name = 'tracer mass' + + def __init__(self, tracer_name, solver_obj, **kwargs): + """ + :arg tracer_name: Name of the tracer. Use canonical field names as in :class:`.FieldDict`. + :arg solver_obj: Thetis solver object + :arg **kwargs: any additional keyword arguments, see DiagnosticCallback + """ + self.name = tracer_name + ' mass' # override name for given tracer + + def mass(): + return comp_tracer_mass_2d(self.solver_obj.fields.elev_2d, + self.solver_obj.fields.bathymetry_2d, + self.solver_obj.fields[tracer_name]) + super(TracerMassConservation2DCallback, self).__init__(mass, solver_obj, **kwargs) + + class TracerMassConservationCallback(ScalarConservationCallback): """Checks conservation of total tracer mass""" name = 'tracer mass' diff --git a/thetis/coupled_timeintegrator_2d.py b/thetis/coupled_timeintegrator_2d.py new file mode 100644 index 000000000..8252397a9 --- /dev/null +++ b/thetis/coupled_timeintegrator_2d.py @@ -0,0 +1,147 @@ +""" +Time integrators for solving coupled shallow water equations with one tracer. +""" +from __future__ import absolute_import +from .utility import * +from . import timeintegrator +from .log import * +from abc import ABCMeta, abstractproperty + + +class CoupledTimeIntegrator2D(timeintegrator.TimeIntegratorBase): + """ + Base class of time integrator for coupled shallow water and tracer equations + """ + __metaclass__ = ABCMeta + + @abstractproperty + def swe_integrator(self): + """time integrator for the shallow water equations""" + pass + + @abstractproperty + def tracer_integrator(self): + """time integrator for the tracer equation""" + pass + + def __init__(self, solver): + """ + :arg solver: :class:`.FlowSolver` object + """ + self.solver = solver + self.options = solver.options + self.fields = solver.fields + self.timesteppers = AttrDict() + print_output('Coupled time integrator: {:}'.format(self.__class__.__name__)) + print_output(' Shallow Water time integrator: {:}'.format(self.swe_integrator.__name__)) + print_output(' Tracer time integrator: {:}'.format(self.tracer_integrator.__name__)) + self._initialized = False + + self._create_integrators() + + def _create_swe_integrator(self): + """ + Create time integrator for 2D system + """ + solver = self.solver + fields = { + 'linear_drag_coefficient': self.options.linear_drag_coefficient, + 'quadratic_drag_coefficient': self.options.quadratic_drag_coefficient, + 'manning_drag_coefficient': self.options.manning_drag_coefficient, + 'viscosity_h': self.options.horizontal_viscosity, + 'lax_friedrichs_velocity_scaling_factor': self.options.lax_friedrichs_velocity_scaling_factor, + 'coriolis': self.options.coriolis_frequency, + 'wind_stress': self.options.wind_stress, + 'atmospheric_pressure': self.options.atmospheric_pressure, + 'momentum_source': self.options.momentum_source_2d, + 'volume_source': self.options.volume_source_2d, } + + if issubclass(self.swe_integrator, timeintegrator.CrankNicolson): + self.timesteppers.swe2d = self.swe_integrator( + solver.eq_sw, self.fields.solution_2d, + fields, solver.dt, + bnd_conditions=solver.bnd_functions['shallow_water'], + solver_parameters=self.options.timestepper_options.solver_parameters, + semi_implicit=self.options.timestepper_options.use_semi_implicit_linearization, + theta=self.options.timestepper_options.implicitness_theta) + else: + self.timesteppers.swe2d = self.swe_integrator(solver.eq_sw, self.fields.solution_2d, + fields, solver.dt, + bnd_conditions=solver.bnd_functions['shallow_water'], + solver_parameters=self.options.timestepper_options.solver_parameters) + + def _create_tracer_integrator(self): + """ + Create time integrator for tracer equation + """ + solver = self.solver + + if self.solver.options.solve_tracer: + uv, elev = self.fields.solution_2d.split() + fields = {'elev_2d': elev, + 'uv_2d': uv, + 'diffusivity_h': self.options.horizontal_diffusivity, + 'source': self.options.tracer_source_2d, + 'lax_friedrichs_tracer_scaling_factor': self.options.lax_friedrichs_tracer_scaling_factor, + } + if issubclass(self.tracer_integrator, timeintegrator.CrankNicolson): + self.timesteppers.tracer = self.tracer_integrator( + solver.eq_tracer, solver.fields.tracer_2d, fields, solver.dt, + bnd_conditions=solver.bnd_functions['tracer'], + solver_parameters=self.options.timestepper_options.solver_parameters_tracer, + semi_implicit=self.options.timestepper_options.use_semi_implicit_linearization, + theta=self.options.timestepper_options.implicitness_theta) + else: + self.timesteppers.tracer = self.tracer_integrator(solver.eq_tracer, solver.fields.tracer_2d, fields, solver.dt, + bnd_conditions=solver.bnd_functions['tracer'], + solver_parameters=self.options.timestepper_options.solver_parameters_tracer,) + + def _create_integrators(self): + """ + Creates all time integrators with the correct arguments + """ + self._create_swe_integrator() + self._create_tracer_integrator() + self.cfl_coeff_2d = min(self.timesteppers.swe2d.cfl_coeff, self.timesteppers.tracer.cfl_coeff) + + def set_dt(self, dt): + """ + Set time step for the coupled time integrator + + :arg float dt: Time step. + """ + for stepper in sorted(self.timesteppers): + self.timesteppers[stepper].set_dt(dt) + + def initialize(self, solution2d): + """ + Assign initial conditions to all necessary fields + + Initial conditions are read from :attr:`fields` dictionary. + """ + + # solution2d is only provided to make a timeintegrator of just the swe + # compatible with the 2d coupled timeintegrator + assert solution2d == self.fields.solution_2d + + self.timesteppers.swe2d.initialize(self.fields.solution_2d) + if self.options.solve_tracer: + self.timesteppers.tracer.initialize(self.fields.tracer_2d) + + self._initialized = True + + def advance(self, t, update_forcings=None): + self.timesteppers.swe2d.advance(t, update_forcings=update_forcings) + self.timesteppers.tracer.advance(t, update_forcings=update_forcings) + if self.options.use_limiter_for_tracers: + self.solver.tracer_limiter.apply(self.fields.tracer_2d) + + +class CoupledCrankNicolson2D(CoupledTimeIntegrator2D): + swe_integrator = timeintegrator.CrankNicolson + tracer_integrator = timeintegrator.CrankNicolson + + +class CoupledCrankEuler2D(CoupledTimeIntegrator2D): + swe_integrator = timeintegrator.CrankNicolson + tracer_integrator = timeintegrator.ForwardEuler diff --git a/thetis/field_defs.py b/thetis/field_defs.py index 6821527fc..7a9c39919 100644 --- a/thetis/field_defs.py +++ b/thetis/field_defs.py @@ -64,6 +64,12 @@ 'unit': 'm s-1', 'filename': 'Velocity2d', } +field_metadata['tracer_2d'] = { + 'name': 'Depth averaged tracer', + 'shortname': 'Tracer', + 'unit': '', + 'filename': 'Tracer2d', +} field_metadata['uv_dav_2d'] = { 'name': 'Depth averaged velocity', 'shortname': 'Depth averaged velocity', diff --git a/thetis/options.py b/thetis/options.py index 8d20f8e5e..95fe1a9ff 100644 --- a/thetis/options.py +++ b/thetis/options.py @@ -25,6 +25,10 @@ class SemiImplicitTimestepperOptions2d(TimeStepperOptions): 'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'multiplicative', }).tag(config=True) + solver_parameters_tracer = PETScSolverParameters({ + 'ksp_type': 'gmres', + 'pc_type': 'sor', + }).tag(config=True) class SteadyStateTimestepperOptions2d(TimeStepperOptions): @@ -351,6 +355,13 @@ class CommonModelOptions(FrozenConfigurable): True, help="use Lax Friedrichs stabilisation in horizontal momentum advection.").tag(config=True) lax_friedrichs_velocity_scaling_factor = FiredrakeConstant( Constant(1.0), help="Scaling factor for Lax Friedrichs stabilisation term in horiozonal momentum advection.").tag(config=True) + use_lax_friedrichs_tracer = Bool( + True, help="Use Lax Friedrichs stabilisation in tracer advection.").tag(config=True) + lax_friedrichs_tracer_scaling_factor = FiredrakeConstant( + Constant(1.0), help="Scaling factor for tracer Lax Friedrichs stability term.").tag(config=True) + use_limiter_for_tracers = Bool( + True, help="Apply P1DG limiter for tracer fields").tag(config=True) + check_volume_conservation_2d = Bool( False, help=""" Compute volume of the 2D mode at every export @@ -437,6 +448,10 @@ class CommonModelOptions(FrozenConfigurable): None, allow_none=True, help="Source term for 2D momentum equation").tag(config=True) volume_source_2d = FiredrakeScalarExpression( None, allow_none=True, help="Source term for 2D continuity equation").tag(config=True) + tracer_source_2d = FiredrakeScalarExpression( + None, allow_none=True, help="Source term for 2D tracer equation").tag(config=True) + horizontal_diffusivity = FiredrakeCoefficient( + None, allow_none=True, help="Horizontal diffusivity for tracers").tag(config=True) # NOTE all parameters are now case sensitive @@ -459,6 +474,7 @@ class CommonModelOptions(FrozenConfigurable): class ModelOptions2d(CommonModelOptions): """Options for 2D depth-averaged shallow water model""" name = 'Depth-averaged 2D model' + solve_tracer = Bool(False, help='Solve tracer transport').tag(config=True) use_wetting_and_drying = Bool( False, help=r"""bool: Turn on wetting and drying @@ -472,6 +488,20 @@ class ModelOptions2d(CommonModelOptions): Used in bathymetry displacement function that ensures positive water depths. Unit is meters. """).tag(config=True) + check_tracer_conservation = Bool( + False, help=""" + Compute total tracer mass at every export + + Prints deviation from the initial mass to stdout. + """).tag(config=True) + + check_tracer_overshoot = Bool( + False, help=""" + Compute tracer overshoots at every export + + Prints overshoot values that exceed the initial range to stdout. + """).tag(config=True) + @attach_paired_options("timestepper_type", PairedEnum([('LeapFrog', ExplicitTimestepperOptions3d), @@ -527,8 +557,6 @@ class ModelOptions3d(CommonModelOptions): See :class:`.SmagorinskyViscosity`.""").tag(config=True) - use_limiter_for_tracers = Bool( - True, help="Apply P1DG limiter for tracer fields").tag(config=True) use_limiter_for_velocity = Bool( True, help="Apply P1DG limiter for 3D horizontal velocity field").tag(config=True) use_lax_friedrichs_tracer = Bool( diff --git a/thetis/solver2d.py b/thetis/solver2d.py index 75fd61eed..bd95e3836 100644 --- a/thetis/solver2d.py +++ b/thetis/solver2d.py @@ -7,6 +7,9 @@ from . import timeintegrator from . import rungekutta from . import implicitexplicit +from . import coupled_timeintegrator_2d +from . import tracer_eq_2d +import weakref import time as time_mod from mpi4py import MPI from . import exporter @@ -15,6 +18,7 @@ from . import callback from .log import * from collections import OrderedDict +import thetis.limiter as limiter class FlowSolver2d(FrozenClass): @@ -124,7 +128,7 @@ def __init__(self, mesh2d, bathymetry_2d, options=None): self.export_initial_state = True """Do export initial state. False if continuing a simulation""" - self.bnd_functions = {'shallow_water': {}} + self.bnd_functions = {'shallow_water': {}, 'tracer': {}} self._isfrozen = True @@ -212,6 +216,8 @@ def create_function_spaces(self): raise Exception('Unsupported finite element family {:}'.format(self.options.element_family)) self.function_spaces.V_2d = MixedFunctionSpace([self.function_spaces.U_2d, self.function_spaces.H_2d]) + self.function_spaces.Q_2d = FunctionSpace(self.mesh2d, 'DG', 1, name='Q_2d') + self._isfrozen = True def create_equations(self): @@ -233,6 +239,15 @@ def create_equations(self): self.options ) self.eq_sw.bnd_functions = self.bnd_functions['shallow_water'] + if self.options.solve_tracer: + self.fields.tracer_2d = Function(self.function_spaces.Q_2d, name='tracer_2d') + self.eq_tracer = tracer_eq_2d.TracerEquation2D(self.function_spaces.Q_2d, bathymetry=self.fields.bathymetry_2d, + use_lax_friedrichs=self.options.use_lax_friedrichs_tracer) + if self.options.use_limiter_for_tracers and self.options.polynomial_degree > 0: + self.tracer_limiter = limiter.VertexBasedP1DGLimiter(self.function_spaces.Q_2d) + else: + self.tracer_limiter = None + self._isfrozen = True # disallow creating new attributes def create_timestepper(self): @@ -279,12 +294,15 @@ def create_timestepper(self): bnd_conditions=self.bnd_functions['shallow_water'], solver_parameters=self.options.timestepper_options.solver_parameters) elif self.options.timestepper_type == 'CrankNicolson': - self.timestepper = timeintegrator.CrankNicolson(self.eq_sw, self.fields.solution_2d, - fields, self.dt, - bnd_conditions=self.bnd_functions['shallow_water'], - solver_parameters=self.options.timestepper_options.solver_parameters, - semi_implicit=self.options.timestepper_options.use_semi_implicit_linearization, - theta=self.options.timestepper_options.implicitness_theta) + if self.options.solve_tracer: + self.timestepper = coupled_timeintegrator_2d.CoupledCrankNicolson2D(weakref.proxy(self)) + else: + self.timestepper = timeintegrator.CrankNicolson(self.eq_sw, self.fields.solution_2d, + fields, self.dt, + bnd_conditions=self.bnd_functions['shallow_water'], + solver_parameters=self.options.timestepper_options.solver_parameters, + semi_implicit=self.options.timestepper_options.use_semi_implicit_linearization, + theta=self.options.timestepper_options.implicitness_theta) elif self.options.timestepper_type == 'DIRK22': self.timestepper = rungekutta.DIRK22(self.eq_sw, self.fields.solution_2d, fields, self.dt, @@ -384,7 +402,7 @@ def initialize(self): self.create_exporters() self._initialized = True - def assign_initial_conditions(self, elev=None, uv=None): + def assign_initial_conditions(self, elev=None, uv=None, tracer=None): """ Assigns initial conditions @@ -400,6 +418,8 @@ def assign_initial_conditions(self, elev=None, uv=None): elev_2d.project(elev) if uv is not None: uv_2d.project(uv) + if tracer is not None and self.options.solve_tracer: + self.fields.tracer_2d.project(tracer) self.timestepper.initialize(self.fields.solution_2d) @@ -515,6 +535,8 @@ def iterate(self, update_forcings=None, if not self._initialized: self.initialize() + self.options.use_limiter_for_tracers &= self.options.polynomial_degree > 0 + t_epsilon = 1.0e-5 cputimestamp = time_mod.clock() next_export_t = self.simulation_time + self.options.simulation_export_time @@ -526,6 +548,20 @@ def iterate(self, update_forcings=None, append_to_log=True) self.add_callback(c) + if self.options.check_tracer_conservation: + c = callback.TracerMassConservation2DCallback('tracer_2d', + self, + export_to_hdf5=dump_hdf5, + append_to_log=True) + self.add_callback(c, eval_interval='export') + + if self.options.check_tracer_overshoot: + c = callback.TracerOvershootCallBack('tracer_2d', + self, + export_to_hdf5=dump_hdf5, + append_to_log=True) + self.add_callback(c, eval_interval='export') + # initial export self.print_state(0.0) if self.export_initial_state: diff --git a/thetis/tracer_eq_2d.py b/thetis/tracer_eq_2d.py new file mode 100644 index 000000000..7df112382 --- /dev/null +++ b/thetis/tracer_eq_2d.py @@ -0,0 +1,270 @@ +r""" +2D advection diffusion equation for tracers. + +The advection-diffusion equation of tracer :math:`T` in non-conservative form reads + +.. math:: + \frac{\partial T}{\partial t} + + \nabla_h \cdot (\textbf{u} T) + = \nabla_h \cdot (\mu_h \nabla_h T) + :label: tracer_eq + +where :math:`\nabla_h` denotes horizontal gradient, :math:`\textbf{u}` are the horizontal +velocities, respectively, and +:math:`\mu_h` denote horizontal and vertical diffusivity. +""" +from __future__ import absolute_import +from .utility import * +from .equation import Term, Equation + +__all__ = [ + 'TracerEquation2D', + 'TracerTerm', + 'HorizontalAdvectionTerm', + 'HorizontalDiffusionTerm', + 'SourceTerm', +] + + +class TracerTerm(Term): + """ + Generic tracer term that provides commonly used members and mapping for + boundary functions. + """ + def __init__(self, function_space, + bathymetry=None, use_lax_friedrichs=True): + """ + :arg function_space: :class:`FunctionSpace` where the solution belongs + :kwarg bathymetry: bathymetry of the domain + :type bathymetry: 2D :class:`Function` or :class:`Constant` + """ + super(TracerTerm, self).__init__(function_space) + self.bathymetry = bathymetry + self.cellsize = CellSize(self.mesh) + continuity = element_continuity(self.function_space.ufl_element()) + self.horizontal_dg = continuity.horizontal == 'dg' + self.use_lax_friedrichs = use_lax_friedrichs + + # define measures with a reasonable quadrature degree + p = self.function_space.ufl_element().degree() + self.quad_degree = 2*p + 1 + self.dx = dx(degree=self.quad_degree) + self.dS = dS(degree=self.quad_degree) + self.ds = ds(degree=self.quad_degree) + + def get_bnd_functions(self, c_in, uv_in, elev_in, bnd_id, bnd_conditions): + """ + Returns external values of tracer and uv for all supported + boundary conditions. + + Volume flux (flux) and normal velocity (un) are defined positive out of + the domain. + + :arg c_in: Internal value of tracer + :arg uv_in: Internal value of horizontal velocity + :arg elev_in: Internal value of elevation + :arg bnd_id: boundary id + :type bnd_id: int + :arg bnd_conditions: dict of boundary conditions: + ``{bnd_id: {field: value, ...}, ...}`` + """ + funcs = bnd_conditions.get(bnd_id) + + if 'elev' in funcs: + elev_ext = funcs['elev'] + else: + elev_ext = elev_in + if 'value' in funcs: + c_ext = funcs['value'] + else: + c_ext = c_in + if 'uv' in funcs: + uv_ext = funcs['uv'] + elif 'flux' in funcs: + assert self.bathymetry is not None + h_ext = elev_ext + self.bathymetry + area = h_ext*self.boundary_len # NOTE using external data only + uv_ext = funcs['flux']/area*self.normal + elif 'un' in funcs: + uv_ext = funcs['un']*self.normal + else: + uv_ext = uv_in + + return c_ext, uv_ext, elev_ext + + +class HorizontalAdvectionTerm(TracerTerm): + r""" + Advection of tracer term, :math:`\bar{\textbf{u}} \cdot \nabla T` + + The weak form is + + .. math:: + \int_\Omega \bar{\textbf{u}} \cdot \boldsymbol{\psi} \cdot \nabla T dx + = - \int_\Omega \nabla_h \cdot (\bar{\textbf{u}} \boldsymbol{\psi}) \cdot T dx + + \int_\Gamma \text{avg}(T) \cdot \text{jump}(\boldsymbol{\psi} + (\bar{\textbf{u}}\cdot\textbf{n})) dS + + where the right hand side has been integrated by parts; + :math:`\textbf{n}` is the unit normal of + the element interfaces, and :math:`\text{jump}` and :math:`\text{avg}` denote the + jump and average operators across the interface. + """ + def residual(self, solution, solution_old, fields, fields_old, bnd_conditions=None): + if fields_old.get('uv_2d') is None: + return 0 + elev = fields_old['elev_2d'] + uv = fields_old['uv_2d'] + + uv_p1 = fields_old.get('uv_p1') + uv_mag = fields_old.get('uv_mag') + # FIXME is this an option? + lax_friedrichs_factor = fields_old.get('lax_friedrichs_tracer_scaling_factor') + + f = 0 + f += -(Dx(uv[0] * self.test, 0) * solution + + Dx(uv[1] * self.test, 1) * solution) * self.dx + + if self.horizontal_dg: + # add interface term + uv_av = avg(uv) + un_av = (uv_av[0]*self.normal('-')[0] + + uv_av[1]*self.normal('-')[1]) + s = 0.5*(sign(un_av) + 1.0) + c_up = solution('-')*s + solution('+')*(1-s) + + f += c_up*(jump(self.test, uv[0] * self.normal[0]) + + jump(self.test, uv[1] * self.normal[1])) * self.dS + # Lax-Friedrichs stabilization + if self.use_lax_friedrichs: + if uv_p1 is not None: + gamma = 0.5*abs((avg(uv_p1)[0]*self.normal('-')[0] + + avg(uv_p1)[1]*self.normal('-')[1]))*lax_friedrichs_factor + elif uv_mag is not None: + gamma = 0.5*avg(uv_mag)*lax_friedrichs_factor + else: + gamma = 0.5*abs(un_av)*lax_friedrichs_factor + f += gamma*dot(jump(self.test), jump(solution))*self.dS + if bnd_conditions is not None: + for bnd_marker in self.boundary_markers: + funcs = bnd_conditions.get(bnd_marker) + ds_bnd = ds(int(bnd_marker), degree=self.quad_degree) + c_in = solution + if funcs is None: + f += c_in * (uv[0]*self.normal[0] + + uv[1]*self.normal[1])*self.test*ds_bnd + else: + c_ext, uv_ext, eta_ext = self.get_bnd_functions(c_in, uv, elev, bnd_marker, bnd_conditions) + uv_av = 0.5*(uv + uv_ext) + un_av = self.normal[0]*uv_av[0] + self.normal[1]*uv_av[1] + s = 0.5*(sign(un_av) + 1.0) + c_up = c_in*s + c_ext*(1-s) + f += c_up*(uv_av[0]*self.normal[0] + + uv_av[1]*self.normal[1])*self.test*ds_bnd + + return -f + + +class HorizontalDiffusionTerm(TracerTerm): + r""" + Horizontal diffusion term :math:`-\nabla_h \cdot (\mu_h \nabla_h T)` + + Using the symmetric interior penalty method the weak form becomes + + .. math:: + -\int_\Omega \nabla_h \cdot (\mu_h \nabla_h T) \phi dx + =& \int_\Omega \mu_h (\nabla_h \phi) \cdot (\nabla_h T) dx \\ + &- \int_{\mathcal{I}_h\cup\mathcal{I}_v} \text{jump}(\phi \textbf{n}_h) + \cdot \text{avg}(\mu_h \nabla_h T) dS + - \int_{\mathcal{I}_h\cup\mathcal{I}_v} \text{jump}(T \textbf{n}_h) + \cdot \text{avg}(\mu_h \nabla \phi) dS \\ + &+ \int_{\mathcal{I}_h\cup\mathcal{I}_v} \sigma \text{avg}(\mu_h) \text{jump}(T \textbf{n}_h) \cdot + \text{jump}(\phi \textbf{n}_h) dS + + where :math:`\sigma` is a penalty parameter, + see Epshteyn and Riviere (2007). + + Epshteyn and Riviere (2007). Estimation of penalty parameters for symmetric + interior penalty Galerkin methods. Journal of Computational and Applied + Mathematics, 206(2):843-872. http://dx.doi.org/10.1016/j.cam.2006.08.029 + + """ + def residual(self, solution, solution_old, fields, fields_old, bnd_conditions=None): + if fields_old.get('diffusivity_h') is None: + return 0 + diffusivity_h = fields_old['diffusivity_h'] + diff_tensor = as_matrix([[diffusivity_h, 0, ], + [0, diffusivity_h, ]]) + grad_test = grad(self.test) + diff_flux = dot(diff_tensor, grad(solution)) + + f = 0 + f += inner(grad_test, diff_flux)*self.dx + + if self.horizontal_dg: + # Interior Penalty method by + # Epshteyn (2007) doi:10.1016/j.cam.2006.08.029 + # sigma = 3*k_max**2/k_min*p*(p+1)*cot(Theta) + # k_max/k_min - max/min diffusivity + # p - polynomial degree + # Theta - min angle of triangles + # assuming k_max/k_min=2, Theta=pi/3 + # sigma = 6.93 = 3.5*p*(p+1) + + degree_h = self.function_space.ufl_element().degree() + sigma = 5.0*degree_h*(degree_h + 1)/self.cellsize + if degree_h == 0: + sigma = 1.5 / self.cellsize + alpha = avg(sigma) + ds_interior = self.dS + f += alpha*inner(jump(self.test, self.normal), + dot(avg(diff_tensor), jump(solution, self.normal)))*ds_interior + f += -inner(avg(dot(diff_tensor, grad(self.test))), + jump(solution, self.normal))*ds_interior + f += -inner(jump(self.test, self.normal), + avg(dot(diff_tensor, grad(solution))))*ds_interior + + return -f + + +class SourceTerm(TracerTerm): + """ + Generic source term + + The weak form reads + + .. math:: + F_s = \int_\Omega \sigma \phi dx + + where :math:`\sigma` is a user defined scalar :class:`Function`. + + """ + def residual(self, solution, solution_old, fields, fields_old, bnd_conditions=None): + f = 0 + source = fields_old.get('source') + if source is not None: + f += -inner(source, self.test)*self.dx + return -f + + +class TracerEquation2D(Equation): + """ + 2D tracer advection-diffusion equation :eq:`tracer_eq` in conservative form + """ + def __init__(self, function_space, + bathymetry=None, + use_lax_friedrichs=False): + """ + :arg function_space: :class:`FunctionSpace` where the solution belongs + :kwarg bathymetry: bathymetry of the domain + :type bathymetry: 2D :class:`Function` or :class:`Constant` + + :kwarg bool use_symmetric_surf_bnd: If True, use symmetric surface boundary + condition in the horizontal advection term + """ + super(TracerEquation2D, self).__init__(function_space) + + args = (function_space, bathymetry, use_lax_friedrichs) + self.add_term(HorizontalAdvectionTerm(*args), 'explicit') + self.add_term(HorizontalDiffusionTerm(*args), 'explicit') + self.add_term(SourceTerm(*args), 'source') diff --git a/thetis/utility.py b/thetis/utility.py index 74cd5e9e6..dd4a8adf5 100644 --- a/thetis/utility.py +++ b/thetis/utility.py @@ -274,6 +274,18 @@ def comp_volume_3d(mesh): return val +def comp_tracer_mass_2d(eta, bath, scalar_func): + """ + Computes total tracer mass in the 2D domain + :arg eta: elevation :class: 'Function' + :arg bath: bathymetry :class: 'Function' + :arg scalar_func: scalar :class:`Function` to integrate + """ + + val = assemble((eta+bath)*scalar_func*dx) + return val + + def comp_tracer_mass_3d(scalar_func): """ Computes total tracer mass in the 3D domain