From 477ef7724b9d623b185bd3e3c2f79dde936455f0 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 1 Nov 2018 17:04:31 -0500 Subject: [PATCH 1/2] Add testing suite for C++ code generation in Dagrt --- requirements.txt | 2 +- test/lapack_kernels.H | 12 + test/sim_types.H | 22 ++ test/test_ab_squarewave.c | 93 ++++++ test/test_codegen_cxx.py | 592 ++++++++++++++++++++++++++++++++++ test/test_fancy_rk.c | 68 ++++ test/test_mrab.c | 121 +++++++ test/test_mrab_squarewave.c | 105 ++++++ test/test_rk.c | 79 +++++ test/test_rk_adaptive.c | 94 ++++++ test/test_rk_adaptive_error.c | 67 ++++ test/timing.H | 18 ++ 12 files changed, 1272 insertions(+), 1 deletion(-) create mode 100644 test/lapack_kernels.H create mode 100644 test/sim_types.H create mode 100644 test/test_ab_squarewave.c create mode 100755 test/test_codegen_cxx.py create mode 100644 test/test_fancy_rk.c create mode 100644 test/test_mrab.c create mode 100644 test/test_mrab_squarewave.c create mode 100644 test/test_rk.c create mode 100644 test/test_rk_adaptive.c create mode 100644 test/test_rk_adaptive_error.c create mode 100644 test/timing.H diff --git a/requirements.txt b/requirements.txt index e22f333..759abbb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ git+git://github.com/inducer/pytools git+git://github.com/inducer/pymbolic -git+https://gitlab.tiker.net/inducer/dagrt.git \ No newline at end of file +git+https://gitlab.tiker.net/inducer/dagrt.git @cory_cxx_generator diff --git a/test/lapack_kernels.H b/test/lapack_kernels.H new file mode 100644 index 0000000..e48dda7 --- /dev/null +++ b/test/lapack_kernels.H @@ -0,0 +1,12 @@ +#ifndef __LAPACK_KERNELS_H__ +#define __LAPACK_KERNELS_H__ + +//#include "FC.h" + +extern "C" { + + void dgesv_(int*, int*, double*, int*, int*, double*, int*, int*); + void dgesvd_(char*, char*, int*, int*, double [], int*, double*, double*, int*, double*, int*, double [], int*, int*); + +} +#endif diff --git a/test/sim_types.H b/test/sim_types.H new file mode 100644 index 0000000..c7a2eb1 --- /dev/null +++ b/test/sim_types.H @@ -0,0 +1,22 @@ +// sim types +#include +#include +#include +#include +#include +#include +#include +using namespace std; + +struct region_type { + + int n_grids; + int* n_grid_dofs; + +}; + +struct sim_grid_state_type { + + std::shared_ptr conserved_var; + +}; diff --git a/test/test_ab_squarewave.c b/test/test_ab_squarewave.c new file mode 100644 index 0000000..7cb5196 --- /dev/null +++ b/test/test_ab_squarewave.c @@ -0,0 +1,93 @@ +#include +#include +#include +#include +#include +#include "lapack_kernels.H" +#include "sim_types.H" +#include "timing.H" +#include "ABMethod.H" +using namespace std; + +int main() +{ + + dagrt_state_type state; + + std::shared_ptr initial_condition(new double [2]); + std::shared_ptr t_ptr(new double); + std::shared_ptr dt_ptr(new double); + double true_sol[2]; + int ntrips[2]; + + int run_count=2; + double t_fin=1.0; + double t = 0.0; + + double dt_values[run_count]; + double errors[run_count]; + + double est_order; + double min_order; + + int stderr=0; + + int istep; + int irun; + int k; + + // start code ---------------------------------------------------------------- + + initial_condition.get()[0] = 2; + initial_condition.get()[1] = 2.3; + + ntrips[0] = 20; + ntrips[1] = 50; + + k = 0; + + for(int irun = 0; irun < run_count; irun++) + { + dt_values[irun] = t_fin/ntrips[irun]; + *dt_ptr = dt_values[irun]; + *t_ptr = t; + + initialize( + state, + initial_condition, + t_ptr, + dt_ptr); + + for (int istep = 0; istep < ntrips[irun]; istep++) + { + if (k == 1) + { + state.dagrt_dt = dt_values[irun]/4; + k = 0; + } + else + { + state.dagrt_dt = dt_values[irun]; + k = 1; + } + run(state); + } + + true_sol[0] = initial_condition.get()[0] * exp(-2*state.ret_time_y); + true_sol[1] = initial_condition.get()[1] * exp(-2*state.ret_time_y); + errors[irun] = sqrt((pow((true_sol[0]-state.ret_state_y.get()[0]),2)) + pow((true_sol[1]-state.ret_state_y.get()[1]),2)); + + shutdown(state); + cout << "done" << endl; + } + + min_order = MIN_ORDER; + est_order = log(errors[1]/errors[0])/log(dt_values[1]/dt_values[0]); + + cout << "estimated order: " << est_order << endl; + if (est_order < min_order) + { + cerr << "ERROR: achieved order too low: " << est_order << " < " << min_order << endl; + } + +} diff --git a/test/test_codegen_cxx.py b/test/test_codegen_cxx.py new file mode 100755 index 0000000..e966a2e --- /dev/null +++ b/test/test_codegen_cxx.py @@ -0,0 +1,592 @@ +#! /usr/bin/env python +from __future__ import division, with_statement, print_function + +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner, Matt Wala" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import sys +import pytest + +import dagrt.codegen.cxx as cxx +from leap.rk import ODE23Method, ODE45Method, RK4Method, LSRK4Method + +from leap.multistep.multirate import TwoRateAdamsBashforthMethod + +from dagrt.utils import run_cxx + +#skip = pytest.mark.skipif(True, reason="not fully implemented") + + +def read_file(rel_path): + from os.path import join, abspath, dirname + path = join(abspath(dirname(__file__)), rel_path) + with open(path, "r") as inf: + return inf.read() + +# {{{ test rk methods + + +@pytest.mark.parametrize(("min_order", "stepper"), [ + (2, ODE23Method("y", use_high_order=False)), + (3, ODE23Method("y", use_high_order=True)), + (4, ODE45Method("y", use_high_order=False)), + (5, ODE45Method("y", use_high_order=True)), + (4, RK4Method("y")), + (4, LSRK4Method("y")), + ]) +def test_rk_codegen(min_order, stepper, print_code=False): + """Test whether CXX code generation for the Runge-Kutta + timestepper works. + """ + + component_id = 'y' + rhs_function = 'y' + + from dagrt.function_registry import ( + base_function_registry, register_ode_rhs) + freg = register_ode_rhs(base_function_registry, component_id, + identifier=rhs_function) + freg = freg.register_codegen(rhs_function, "cxx", + cxx.CallCode(""" + for (int i = 0; i < 2; i++) + { + ${result}[i] = -2 * ${y}[i]; + } + """)) + + code = stepper.generate() + + codegen = cxx.CodeGenerator( + 'RKMethod', + user_type_map={ + component_id: cxx.ArrayType( + (2,), + cxx.BuiltinType('double'), + ) + }, + function_registry=freg, + ) + + code_str = codegen(code) + if print_code: + print(code_str) + + run_cxx([ + ("test_rk.c", read_file("test_rk.c").replace( + "MIN_ORDER", str(min_order - 0.3))), + ("RKMethod.H", code_str), + ]) + +# }}} + + +# {{{ test fancy codegen + +def test_rk_codegen_fancy(): + """Test whether CXX code generation with lots of fancy features for the + Runge-Kutta timestepper works. + """ + + component_id = 'y' + rhs_function = 'y' + state_filter_name = 'state_filter_y' + + stepper = ODE23Method(component_id, use_high_order=True, + state_filter_name=state_filter_name) + + from dagrt.function_registry import ( + base_function_registry, register_ode_rhs, + register_function, UserType) + freg = register_ode_rhs(base_function_registry, component_id, + identifier=rhs_function) + freg = freg.register_codegen(rhs_function, "cxx", + cxx.CallCode(""" + <% + + igrid = declare_new("int", "igrid") + i = declare_new("int", "i") + + %> + + for (int ${igrid} = 0; ${igrid} < region.n_grids; ${igrid}++) + { + for (int ${i} = 0; ${i} < region.n_grid_dofs[${igrid}]; ${i}++) + { + ${result}[${igrid}].conserved_var.get()[${i}] = + -2 * ${y}[${igrid}].conserved_var.get()[${i}]; + } + } + + """)) + freg = register_function(freg, "notify_pre_state_update", ("updated_component",)) + freg = freg.register_codegen("notify_pre_state_update", "cxx", + cxx.CallCode(""" + cout << "before state update" << endl; + """)) + freg = register_function( + freg, "notify_post_state_update", ("updated_component",)) + freg = freg.register_codegen("notify_post_state_update", "cxx", + cxx.CallCode(""" + cout << "after state update" << endl; + """)) + + freg = register_function(freg, ""+state_filter_name, ("y",), + result_names=("result",), result_kinds=(UserType("y"),)) + freg = freg.register_codegen(""+state_filter_name, "cxx", + cxx.CallCode(""" + // mess with state + <% + + igrid = declare_new("int", "igrid") + i = declare_new("int", "i") + + %> + + for (int ${igrid} = 0; ${igrid} < region.n_grids; ${igrid}++) + { + for (int ${i} = 0; ${i} < region.n_grid_dofs[${igrid}]; ${i}++) + { + ${result}[${igrid}].conserved_var.get()[${i}] = + 0.95 * ${y}[${igrid}].conserved_var.get()[${i}]; + } + } + + """)) + + code = stepper.generate() + + codegen = cxx.CodeGenerator( + 'RKMethod', + user_type_map={ + component_id: cxx.ArrayType( + "region.n_grids", + index_vars="igrid", + element_type=cxx.StructureType( + "sim_grid_state_type", + ( + ("conserved_var", cxx.PointerType( + cxx.ArrayType( + ("region.n_grid_dofs[igrid]",), + cxx.BuiltinType('double')))), + ))) + }, + function_registry=freg, + call_before_state_update="notify_pre_state_update", + call_after_state_update="notify_post_state_update", + extra_arguments="region", + extra_argument_types="region_type", + emit_instrumentation=True, + timing_function="clock") + + code_str = codegen(code) + with open("RKMethod.H", "w") as outf: + print(code_str, file=outf) + + run_cxx([ + ("test_fancy_rk.c", read_file("test_fancy_rk.c")), + ("RKMethod.H", read_file("RKMethod.H")), + ("sim_types.H", read_file("sim_types.H")), + ("lapack_kernels.H", read_file("lapack_kernels.H")), + ("timing.H", read_file("timing.H")), + ]) + +# }}} + + +@pytest.mark.parametrize("min_order", [2, 3, 4, 5]) +@pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) +def test_multirate_codegen(min_order, method_name): + from leap.multistep.multirate import TwoRateAdamsBashforthMethod + + stepper = TwoRateAdamsBashforthMethod( + method_name, min_order, 4, + # should pass with either, let's alternate by order + # static_dt=True is 'more finnicky', so use that at min_order=5. + static_dt=True if min_order % 2 == 1 else False, + hist_consistency_threshold=1e-8, + early_hist_consistency_threshold="1.25e3 *
**%d" % min_order) + + # Early consistency threshold checked for convergence + # with timestep change - C. Mikida, 2/6/18 (commit hash 2e6ca077) + + # With method 5-Fqs (limiting case), the following maximum relative + # errors were observed: + # for dt = 0.0384: 1.03E-04 + # for dt = 0.0128: 5.43E-08 + + # Reported relative errors motivate constant factor of 1.25e3 for early + # consistency threshold + + code = stepper.generate() + + from dagrt.function_registry import ( + base_function_registry, register_ode_rhs) + + freg = base_function_registry + for func_name in [ + "s2s", + "f2s", + "s2f", + "f2f", + ]: + component_id = { + "s": "slow", + "f": "fast", + }[func_name[-1]] + freg = register_ode_rhs(freg, identifier=func_name, + output_type_id=component_id, + input_type_ids=("slow", "fast"), + input_names=("s", "f")) + + freg = freg.register_codegen("s2f", "cxx", + cxx.CallCode(""" + for (int i = 0; i < 1; i++) + { + ${result}[i] = (sin(2.0 * ${t}) - 1.0) * ${s}[i]; + } + """)) + freg = freg.register_codegen("f2s", "cxx", + cxx.CallCode(""" + for (int i = 0; i < 1; i++) + { + ${result}[i] = (sin(2.0 * ${t}) + 1.0) * ${f}[i]; + } + """)) + freg = freg.register_codegen("f2f", "cxx", + cxx.CallCode(""" + for (int i = 0; i < 1; i++) + { + ${result}[i] = cos(2.0 * ${t}) * ${f}[i]; + } + """)) + freg = freg.register_codegen("s2s", "cxx", + cxx.CallCode(""" + for (int i = 0; i < 1; i++) + { + ${result}[i] = -cos(2.0 * ${t}) * ${s}[i]; + } + """)) + + codegen = cxx.CodeGenerator( + 'MRAB', + user_type_map={ + "slow": cxx.ArrayType( + (1,), + cxx.BuiltinType('double'), + ), + "fast": cxx.ArrayType( + (1,), + cxx.BuiltinType('double'), + ) + }, + function_registry=freg, + ) + + code_str = codegen(code) + + with open("MRABMethod.H", "wt") as outf: + outf.write(code_str) + + fac = 130 + if min_order == 5 and method_name in ["Srs", "Ss"]: + pytest.xfail("Srs,Ss do not achieve fifth order convergence") + + num_trips_one = 10*fac + num_trips_two = 30*fac + + run_cxx([ + ("test_mrab.c", ( + read_file("test_mrab.c") + .replace("MIN_ORDER", str(min_order - 0.3)) + .replace("NUM_TRIPS_ONE", str(num_trips_one)) + .replace("NUM_TRIPS_TWO", str(num_trips_two)))), + ("MRABMethod.H", read_file("MRABMethod.H")), + ("sim_types.H", read_file("sim_types.H")), + ("lapack_kernels.H", read_file("lapack_kernels.H")), + ("timing.H", read_file("timing.H")), + ], + cxx_options=["-llapack", "-lblas"]) + + +def test_adaptive_rk_codegen(): + """Test whether CXX code generation for the Runge-Kutta + timestepper works. + """ + + component_id = 'y' + rhs_function = 'y' + + stepper = ODE45Method(component_id, use_high_order=False, rtol=1e-6) + + from dagrt.function_registry import ( + base_function_registry, register_ode_rhs) + freg = register_ode_rhs(base_function_registry, component_id, + identifier=rhs_function) + freg = freg.register_codegen(rhs_function, "cxx", + cxx.CallCode(""" + ${result}[0] = ${y}[1]; + ${result}[1] = -30 * (pow(${y}[0],2) - 1) * ${y}[1] - ${y}[0]; + """)) + + code = stepper.generate() + + codegen = cxx.CodeGenerator( + 'RKMethod', + user_type_map={ + "y": cxx.ArrayType( + (2,), + cxx.BuiltinType('double'), + ), + }, + function_registry=freg) + + code_str = codegen(code) + with open("RKMethod.H", "wt") as outf: + outf.write(code_str) + + run_cxx([ + ("test_rk_adaptive.c", read_file("test_rk_adaptive.c")), + ("RKMethod.H", read_file("RKMethod.H")), + ("sim_types.H", read_file("sim_types.H")), + ("lapack_kernels.H", read_file("lapack_kernels.H")), + ("timing.H", read_file("timing.H")), + ]) + + +def test_adaptive_rk_codegen_error(): + """Test whether CXX code generation for the Runge-Kutta + timestepper works. + """ + + component_id = 'y' + rhs_function = 'y' + + stepper = ODE45Method(component_id, use_high_order=False, atol=1e-6) + + from dagrt.function_registry import ( + base_function_registry, register_ode_rhs) + freg = register_ode_rhs(base_function_registry, component_id, + identifier=rhs_function) + freg = freg.register_codegen(rhs_function, "cxx", + cxx.CallCode(""" + for (int i = 0; i < 2; i++) + { + ${result}[i] = -2 * ${y}[i]; + } + """)) + + code = stepper.generate() + + codegen = cxx.CodeGenerator( + 'RKMethod', + user_type_map={ + component_id: cxx.ArrayType( + (2,), + cxx.BuiltinType('double'), + ) + }, + function_registry=freg) + + code_str = codegen(code) + with open("RKMethod.H", "wt") as outf: + outf.write(code_str) + + run_cxx([ + ("test_rk_adaptive_error.c", read_file("test_rk_adaptive_error.c")), + ("RKMethod.H", read_file("RKMethod.H")), + ("sim_types.H", read_file("sim_types.H")), + ("lapack_kernels.H", read_file("lapack_kernels.H")), + ("timing.H", read_file("timing.H")), + ]) + + +@pytest.mark.parametrize(("min_order", "hist_length"), [(5, 5), (4, 4), (4, 5), + (3, 3), (3, 4), (2, 2), ]) +def test_singlerate_squarewave(min_order, hist_length): + from leap.multistep import AdamsBashforthMethod + + component_id = 'y' + rhs_function = 'y' + + stepper = AdamsBashforthMethod("y", min_order, hist_length=hist_length) + + from dagrt.function_registry import ( + base_function_registry, register_ode_rhs) + freg = register_ode_rhs(base_function_registry, component_id, + identifier=rhs_function) + freg = freg.register_codegen(rhs_function, "cxx", + cxx.CallCode(""" + for (int i = 0; i < 2; i++) + { + ${result}[i] = -2 * ${y}[i]; + } + """)) + + code = stepper.generate() + + codegen = cxx.CodeGenerator( + 'ABMethod', + user_type_map={ + component_id: cxx.ArrayType( + (2,), + cxx.BuiltinType('double'), + ) + }, + function_registry=freg, + ) + + code_str = codegen(code) + with open("ABMethod.H", "wt") as outf: + outf.write(code_str) + + run_cxx([ + ("test_ab_squarewave.c", read_file("test_ab_squarewave.c").replace( + "MIN_ORDER", str(min_order - 0.3))), + ("ABMethod.H", read_file("ABMethod.H")), + ("sim_types.H", read_file("sim_types.H")), + ("lapack_kernels.H", read_file("lapack_kernels.H")), + ("timing.H", read_file("timing.H")), + ], + cxx_options=["-llapack", "-lblas"]) + + +@pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) +@pytest.mark.parametrize(("min_order", "hist_length"), [(4, 4), (4, 5), + (3, 3), (3, 4), (2, 2), ]) +def test_multirate_squarewave(min_order, hist_length, method_name): + stepper = TwoRateAdamsBashforthMethod(method_name, min_order, 4, + hist_consistency_threshold=1e-8, + early_hist_consistency_threshold="3.5e3 *
** %d" % min_order, + hist_length_slow=hist_length, hist_length_fast=hist_length) + + # Early consistency threshold checked for convergence + # with timestep change - C. Mikida, 2/6/18 (commit hash 0fb6148) + + # With method 3-4-Ss (limiting case), the following maximum relative + # errors were observed: + # for dt = 0.0011: 6.96E-08 + # for dt = 0.00073: 5.90E-09 + + # Reported relative errors motivate constant factor of 3.5e3 for early + # consistency threshold + + code = stepper.generate() + + from dagrt.function_registry import ( + base_function_registry, register_ode_rhs) + + freg = base_function_registry + for func_name in [ + "s2s", + "f2s", + "s2f", + "f2f", + ]: + component_id = { + "s": "slow", + "f": "fast", + }[func_name[-1]] + freg = register_ode_rhs(freg, identifier=func_name, + output_type_id=component_id, + input_type_ids=("slow", "fast"), + input_names=("s", "f")) + + freg = freg.register_codegen("s2f", "cxx", + cxx.CallCode(""" + for (int i = 0; i < 1; i++) + { + ${result}[i] = (sin(2 * ${t}) - 1) * ${s}[i]; + } + """)) + freg = freg.register_codegen("f2s", "cxx", + cxx.CallCode(""" + for (int i = 0; i < 1; i++) + { + ${result}[i] = (sin(2 * ${t}) + 1) * ${f}[i]; + } + """)) + freg = freg.register_codegen("f2f", "cxx", + cxx.CallCode(""" + for (int i = 0; i < 1; i++) + { + ${result}[i] = cos(2 * ${t}) * ${f}[i]; + } + """)) + freg = freg.register_codegen("s2s", "cxx", + cxx.CallCode(""" + for (int i = 0; i < 1; i++) + { + ${result}[i] = -cos(2 * ${t}) * ${s}[i]; + } + """)) + + codegen = cxx.CodeGenerator( + 'MRAB', + user_type_map={ + "slow": cxx.ArrayType( + (1,), + cxx.BuiltinType('double'), + ), + "fast": cxx.ArrayType( + (1,), + cxx.BuiltinType('double'), + ) + }, + function_registry=freg, + ) + + code_str = codegen(code) + with open("MRABMethod.H", "wt") as outf: + outf.write(code_str) + + # Build in conditionals to alter the timestep based on order such that all + # tests pass + + if min_order == 2: + fac = 200 + elif min_order == 3: + fac = 100 + else: + fac = 12 + num_trips_one = 100*fac + num_trips_two = 150*fac + + run_cxx([ + ("test_mrab_squarewave.c", ( + read_file("test_mrab_squarewave.c") + .replace("MIN_ORDER", str(min_order - 0.3)) + .replace("NUM_TRIPS_ONE", str(num_trips_one)) + .replace("NUM_TRIPS_TWO", str(num_trips_two)))), + ("MRABMethod.H", read_file("MRABMethod.H")), + ("sim_types.H", read_file("sim_types.H")), + ("lapack_kernels.H", read_file("lapack_kernels.H")), + ("timing.H", read_file("timing.H")), + ], + cxx_options=["-llapack", "-lblas"]) + + +if __name__ == "__main__": + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) diff --git a/test/test_fancy_rk.c b/test/test_fancy_rk.c new file mode 100644 index 0000000..09446d2 --- /dev/null +++ b/test/test_fancy_rk.c @@ -0,0 +1,68 @@ +#include +#include +#include +#include +#include +#include "lapack_kernels.H" +#include "sim_types.H" +#include "timing.H" +#include "RKMethod.H" +using namespace std; + +int main() +{ + + region_type region; + + std::shared_ptr initial_condition; + + dagrt_state_type dagrt_state; + std::shared_ptr t_ptr(new double); + std::shared_ptr dt_ptr(new double); + + double t_fin=1.0; + double t = 0.0; + double dt = t_fin/20; + int ntrips = 20; + int igrid; + int idof; + int istep; + + // start code ---------------------------------------------------------------- + + region.n_grids = 2; + + region.n_grid_dofs = new int [region.n_grids]; + initial_condition.reset(new sim_grid_state_type [region.n_grids]); + + for (int igrid = 0; igrid < region.n_grids; igrid++) + { + region.n_grid_dofs[igrid] = 2; + initial_condition.get()[igrid].conserved_var.reset(new double [region.n_grid_dofs[igrid]]); + + for (int idof = 0; idof < region.n_grid_dofs[igrid]; idof++) + { + initial_condition.get()[igrid].conserved_var.get()[idof] = 1; + } + } + + *dt_ptr = dt; + *t_ptr = t; + + initialize( + region, + dagrt_state, + initial_condition, + t_ptr, + dt_ptr); + + for (int istep = 0; istep < ntrips; istep++) + { + run(region, dagrt_state); + } + + print_profile(dagrt_state); + shutdown(region, dagrt_state); + +} + diff --git a/test/test_mrab.c b/test/test_mrab.c new file mode 100644 index 0000000..5f46d09 --- /dev/null +++ b/test/test_mrab.c @@ -0,0 +1,121 @@ +#include +#include +#include +#include +#include +#include "lapack_kernels.H" +#include "sim_types.H" +#include "timing.H" +#include "MRABMethod.H" +using namespace std; + +int main() +{ + + dagrt_state_type state; + + std::shared_ptr initial_condition_fast(new double); + std::shared_ptr initial_condition_slow(new double); + std::shared_ptr t_ptr(new double); + std::shared_ptr dt_ptr(new double); + double true_sol_fast; + double true_sol_slow; + int ntrips[2]; + + int run_count = 2; + double t_fin = 50; + + double t = 0.00; + double dt_values[run_count]; + double error_slow[run_count]; + double error_fast[run_count]; + + double min_order; + double est_order_fast; + double est_order_slow; + double err_slow; + double err_fast; + double max_err_slow; + double max_err_fast; + double max_val_slow; + double max_val_fast; + + int stderr=0; + + int irun; + + // start code ---------------------------------------------------------------- + + *initial_condition_fast = (exp(0.0)) * cos(0.0); // fast + *initial_condition_slow = (exp(0.0)) * sin(0.0); // slow + + ntrips[0] = NUM_TRIPS_ONE; + ntrips[1] = NUM_TRIPS_TWO; + + for (int irun = 0; irun < run_count; irun++) + { + dt_values[irun] = t_fin/ntrips[irun]; + *dt_ptr = dt_values[irun]; + *t_ptr = t; + + initialize( + state, + initial_condition_fast, + initial_condition_slow, + t_ptr, + dt_ptr); + + max_err_fast = 0; + max_val_fast = 0; + max_err_slow = 0; + max_val_slow = 0; + + do { + run(state); + + true_sol_fast = (exp(state.ret_time_fast)) * cos(state.ret_time_fast); + true_sol_slow = (exp(state.ret_time_slow)) * sin(state.ret_time_slow); + + err_slow = abs(true_sol_slow - state.ret_state_slow.get()[0]); + err_fast = abs(true_sol_fast - state.ret_state_fast.get()[0]); + + max_err_slow = max(max_err_slow, err_slow); + max_err_fast = max(max_err_fast, err_fast); + + max_val_slow = max(max_val_slow, abs(true_sol_slow)); + max_val_fast = max(max_val_fast, abs(true_sol_fast)); + + + } while (state.ret_time_fast <= t_fin - 1e-10); + + cout << "end time slow: " << state.ret_time_slow << endl; + cout << "end time fast: " << state.ret_time_fast << endl; + + error_slow[irun] = max_err_slow/max_val_slow; + error_fast[irun] = max_err_fast/max_val_fast; + + shutdown(state); + + cout << "done " << dt_values[irun] << " " << + "err_slow: " << error_slow[irun] << " " << + "err_fast: " << error_fast[irun] << endl; + + } + + min_order = MIN_ORDER; + est_order_slow = log(error_slow[1]/error_slow[0])/log(dt_values[1]/dt_values[0]); + est_order_fast = log(error_fast[1]/error_fast[0])/log(dt_values[1]/dt_values[0]); + + cout << "estimated order slow: " << est_order_slow << endl; + if (est_order_slow < min_order) + { + cerr << "ERROR: achieved order too low: " << est_order_slow << " < " << min_order << endl; + } + + cout << "estimated order fast: " << est_order_fast << endl; + if (est_order_fast < min_order) + { + cerr << "ERROR: achieved order too low: " << est_order_fast << " < " << min_order << endl; + } + +} diff --git a/test/test_mrab_squarewave.c b/test/test_mrab_squarewave.c new file mode 100644 index 0000000..3a10e5b --- /dev/null +++ b/test/test_mrab_squarewave.c @@ -0,0 +1,105 @@ +#include +#include +#include +#include +#include +#include "lapack_kernels.H" +#include "sim_types.H" +#include "timing.H" +#include "MRABMethod.H" +using namespace std; + +int main() +{ + + dagrt_state_type state; + + std::shared_ptr initial_condition_fast(new double); + std::shared_ptr initial_condition_slow(new double); + std::shared_ptr t_ptr(new double); + std::shared_ptr dt_ptr(new double); + double true_sol_fast; + double true_sol_slow; + int ntrips[2]; + + int run_count=2; + int k; + double t_fin=11.0; + double t = 0.0; + + double dt_values[run_count]; + double error_slow[run_count]; + double error_fast[run_count]; + + double min_order; + double est_order_fast; + double est_order_slow; + + int stderr=0; + + int irun; + + // start code ---------------------------------------------------------------- + + *initial_condition_fast = (exp(0.0))*cos(0.0); // fast + *initial_condition_slow = (exp(0.0))*sin(0.0); // slow + + ntrips[0] = NUM_TRIPS_ONE; + ntrips[1] = NUM_TRIPS_TWO; + + for (int irun = 0; irun < run_count; irun++) + { + dt_values[irun] = t_fin/ntrips[irun]; + *dt_ptr = dt_values[irun]; + *t_ptr = t; + + initialize( + state, + initial_condition_fast, + initial_condition_slow, + t_ptr, + dt_ptr); + + k = 0; + + do { + if (k == 1) + { + state.dagrt_dt = dt_values[irun]/4; + k = 0; + } + else + { + state.dagrt_dt = dt_values[irun]; + k = 1; + } + run(state); + } while (state.ret_time_fast <= t_fin); + + true_sol_fast = (exp(state.ret_time_fast))*cos(state.ret_time_fast); + true_sol_slow = (exp(state.ret_time_slow))*sin(state.ret_time_slow); + + error_slow[irun] = sqrt((pow((true_sol_slow-state.ret_state_slow.get()[0]),2))); + error_fast[irun] = sqrt((pow((true_sol_fast-state.ret_state_fast.get()[0]),2))); + + shutdown(state); + cout << "done " << dt_values[irun] << " " << error_slow[irun] << " " << error_fast[irun] << endl; + } + + min_order = MIN_ORDER; + est_order_slow = log(error_slow[1]/error_slow[0])/log(dt_values[1]/dt_values[0]); + est_order_fast = log(error_fast[1]/error_fast[0])/log(dt_values[1]/dt_values[0]); + + cout << "estimated order slow: " << est_order_slow << endl; + if (est_order_slow < min_order) + { + cerr << "ERROR: achieved order too low: " << est_order_slow << " < " << min_order << endl; + } + + cout << "estimated order fast: " << est_order_fast << endl; + if (est_order_fast < min_order) + { + cerr << "ERROR: achieved order too low: " << est_order_fast << " < " << min_order << endl; + } + +} diff --git a/test/test_rk.c b/test/test_rk.c new file mode 100644 index 0000000..1cce318 --- /dev/null +++ b/test/test_rk.c @@ -0,0 +1,79 @@ +#include +#include +#include +#include +#include +#include "RKMethod.H" +using namespace std; + +int main() +{ + dagrt_state_type state; + + std::shared_ptr initial_condition(new double [2]); + std::shared_ptr t_ptr(new double); + std::shared_ptr dt_ptr(new double); + double true_sol[2]; + int ntrips[2]; + + int run_count = 2; + double t_fin = 1.0; + + double dt_values[run_count]; + double errors[run_count]; + + double t = 0.00; + double est_order; + double min_order; + + int stderr = 0; + + int istep; + int irun; + + // start code ---------------------------------------------------------------- + + initial_condition.get()[0] = 2; + initial_condition.get()[1] = 2.3; + + ntrips[0] = 20; + ntrips[1] = 50; + + for(int irun = 0; irun < run_count; irun++) + { + dt_values[irun] = t_fin/ntrips[irun]; + + *dt_ptr = dt_values[irun]; + *t_ptr = t; + + initialize( + state, + initial_condition, + t_ptr, + dt_ptr); + + for (int istep = 0; istep < ntrips[irun]; istep++) + { + run(state); + } + + true_sol[0] = initial_condition.get()[0] * exp(-2*state.ret_time_y); + true_sol[1] = initial_condition.get()[1] * exp(-2*state.ret_time_y); + errors[irun] = sqrt((pow((true_sol[0]-state.ret_state_y.get()[0]),2)) + pow((true_sol[1]-state.ret_state_y.get()[1]),2)); + + cout << errors[irun] << endl; + + shutdown(state); + cout << "done" << endl; + } + + min_order = MIN_ORDER; + est_order = log(errors[1]/errors[0])/log(dt_values[1]/dt_values[0]); + + cout << "estimated order: " << est_order << endl; + if (est_order < min_order) + { + cerr << "ERROR: achieved order too low: " << est_order << " < " << min_order << endl; + } +} + diff --git a/test/test_rk_adaptive.c b/test/test_rk_adaptive.c new file mode 100644 index 0000000..1627086 --- /dev/null +++ b/test/test_rk_adaptive.c @@ -0,0 +1,94 @@ +#include +#include +#include +#include +#include +#include "lapack_kernels.H" +#include "sim_types.H" +#include "timing.H" +#include "RKMethod.H" +using namespace std; + +int main() +{ + + dagrt_state_type state; + + std::shared_ptr initial_condition(new double [2]); + std::shared_ptr t_ptr(new double); + std::shared_ptr dt_ptr(new double); + double step_sizes[100]; + + double t_fin=1.0; + double t = 0.0; + double dt_value; + double small_step_frac; + double big_step_frac; + double old_time; + + int stderr=0; + + int istep; + int nruns = 100; + int num_big_steps; + int num_small_steps; + + // start code ---------------------------------------------------------------- + + initial_condition.get()[0] = 2; + initial_condition.get()[1] = 0; + + dt_value = 1e-5; + old_time = 0.0; + num_big_steps = 0; + num_small_steps = 0; + *dt_ptr = dt_value; + *t_ptr = t; + + initialize( + state, + initial_condition, + t_ptr, + dt_ptr); + + for (int istep = 0; istep < nruns; istep++) + { + run(state); + + cout << state.ret_time_y << endl; + step_sizes[istep] = state.ret_time_y - old_time; + old_time = state.ret_time_y; + + if (step_sizes[istep] < 0.01) + { + num_small_steps = num_small_steps + 1; + } + else if (step_sizes[istep] > 0.03) + { + num_big_steps = num_big_steps + 1; + } + + } + + shutdown(state); + cout << "done" << endl; + + cout << num_big_steps << endl; + cout << num_small_steps << endl; + cout << nruns << endl; + big_step_frac = double (num_big_steps) / double (nruns); + small_step_frac = double (num_small_steps) / double (nruns); + + if (big_step_frac>=0.16 && small_step_frac<=0.35) + { + cout << "Test passes: big_step_frac = " << big_step_frac << endl; + cout << "Test passes: small_step_frac = " << small_step_frac << endl; + } + else + { + cerr << "Test fails: big_step_frac = " << big_step_frac << endl; + cerr << "Test fails: small_step_frac = " << small_step_frac << endl; + } + +} + diff --git a/test/test_rk_adaptive_error.c b/test/test_rk_adaptive_error.c new file mode 100644 index 0000000..6305183 --- /dev/null +++ b/test/test_rk_adaptive_error.c @@ -0,0 +1,67 @@ +#include +#include +#include +#include +#include +#include "lapack_kernels.H" +#include "sim_types.H" +#include "timing.H" +#include "RKMethod.H" +using namespace std; + +int main() +{ + + dagrt_state_type state; + + std::shared_ptr initial_condition(new double [2]); + std::shared_ptr t_ptr(new double); + std::shared_ptr dt_ptr(new double); + double true_sol[2]; + + double t_fin = 1.0; + double t = 0.0; + int ntrips = 100; + + double dt_value; + double error; + double atol; + + int stderr = 0; + + int istep; + + // start code ---------------------------------------------------------------- + + initial_condition.get()[0] = 2; + initial_condition.get()[1] = 2.3; + + dt_value = 1e-5; + atol = 1e-6; + *dt_ptr = dt_value; + *t_ptr = t; + + initialize( + state, + initial_condition, + t_ptr, + dt_ptr); + + for (int istep = 1; istep < ntrips; istep++) + { + run(state); + } + + true_sol[0] = initial_condition.get()[0] * exp(-2*state.ret_time_y); + true_sol[1] = initial_condition.get()[1] * exp(-2*state.ret_time_y); + error = sqrt((pow((true_sol[0]-state.ret_state_y.get()[0]),2)) + pow((true_sol[1]-state.ret_state_y.get()[1]),2)); + + shutdown(state); + cout << "done" << endl; + + if (error > atol) + { + cerr << "ERROR: Max error exceeds tolerance: " << error << " < " << atol << endl; + } + +} diff --git a/test/timing.H b/test/timing.H new file mode 100644 index 0000000..d35a289 --- /dev/null +++ b/test/timing.H @@ -0,0 +1,18 @@ +// timing +#include +#include +#include +#include +#include +#include +#include +#include +#include +using namespace std; + +void get_time() +{ + using std::chrono::system_clock; + + system_clock::time_point today = system_clock::now(); +} -- GitLab From 707f26cb15e7aba193f09f2e18bdff80c00c86ec Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 1 Nov 2018 20:43:39 -0500 Subject: [PATCH 2/2] Fixed typo in requirements.txt for dagrt linking --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 759abbb..465e855 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ git+git://github.com/inducer/pytools git+git://github.com/inducer/pymbolic -git+https://gitlab.tiker.net/inducer/dagrt.git @cory_cxx_generator +git+https://gitlab.tiker.net/inducer/dagrt.git@cory_cxx_generator -- GitLab