#882

Committed 31 Aug 2023 07:44PM UTC coverage: 41.798% (-44.7%) from 86.546%

Build # #882

Build Type

push

Committed by

Commit Message

Run Details

19442 of 46514 relevant lines covered (41.8%)

126375.38 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.51

/examples/quickstart/1d_wave_equation.cpp

////////////////////////////////////////////////////////////////////////////////
//  Copyright (c)      2012 Zach Byerly
//  Copyright (c) 2011-2012 Bryce Adelstein-Lelbach
//
//  SPDX-License-Identifier: BSL-1.0
//  Distributed under the Boost Software License, Version 1.0. (See accompanying
//  file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
////////////////////////////////////////////////////////////////////////////////

//
// This is a program written to evolve in time the equation:
//
// D^2 U / Dt^2 = c^2  D^2 U / Dx^2
//
// The parameter alpha = c*dt/dx must be less than 1 to ensure the stability
//     of the algorithm.
// Discretizing the equation and solving for U(t+dt,x) yields
// alpha^2 * (U(t,x+dx)+U(t,x-dx))+2(1-alpha^2)*U(t,x) - U(t-dt,x)
//
// For the first timestep, we approximate U(t-dt,x) by u(t+dt,x)-2*dt*du/dt(t,x)
//

// Include statements.
#include <hpx/config.hpp>
#if !defined(HPX_COMPUTE_DEVICE_CODE)
#include <hpx/chrono.hpp>
#include <hpx/hpx_init.hpp>
#include <hpx/include/async.hpp>
#include <hpx/iostream.hpp>
#include <hpx/modules/actions_base.hpp>
#include <hpx/modules/async_combinators.hpp>
#include <hpx/modules/format.hpp>

#include <cstddef>
#include <cstdint>
#include <fstream>
#include <iostream>
#include <mutex>
#include <vector>

using hpx::program_options::options_description;
using hpx::program_options::value;
using hpx::program_options::variables_map;

using hpx::id_type;
using hpx::invalid_id;

using hpx::async;
using hpx::future;

using hpx::chrono::high_resolution_timer;

using hpx::find_here;

using hpx::cout;

///////////////////////////////////////////////////////////////////////////////
// Globals.

//double const alpha_squared = 0.25;
double alpha_squared = 0;

// Initialized in hpx_main.
id_type here = invalid_id;
double pi = 0.;
double c = 0.;
double dt = 0.;
double dx = 0.;

// Command line argument.
std::uint64_t nt = 0;
std::uint64_t nx = 0;

struct data
{
    // Default constructor: data d1;
    data()
      : mtx()
      , u_value(0.0)
      , computed(false)
    {
    }

    // Copy constructor: data d1; data d2(d1);
    // We can't copy the mutex, because mutexes are noncopyable.
    data(data const& other)
      : mtx()
      , u_value(other.u_value)
      , computed(other.computed)
    {
    }

    data& operator=(data const& other)
    {
        u_value = other.u_value;
        computed = other.computed;
        return *this;
    }

    hpx::mutex mtx;
    double u_value;
    bool computed;
};

std::vector<std::vector<data>> u;

///////////////////////////////////////////////////////////////////////////////
// Forward declaration of the wave function.
double wave(std::uint64_t t, std::uint64_t x);

// Any global function needs to be wrapped into a plain_action if it should be
// invoked as a HPX-thread.
// This generates the required boilerplate we need for remote invocation.
HPX_PLAIN_ACTION(wave)

double calculate_u_tplus_x(
    double u_t_xplus, double u_t_x, double u_t_xminus, double u_tminus_x)
{
    double u_tplus_x = alpha_squared * (u_t_xplus + u_t_xminus) +
        2.0 * (1 - alpha_squared) * u_t_x - u_tminus_x;
    return u_tplus_x;
}

double calculate_u_tplus_x_1st(
    double u_t_xplus, double u_t_x, double u_t_xminus, double u_dot)
{
    double u_tplus_x = 0.5 * alpha_squared * (u_t_xplus + u_t_xminus) +
        (1 - alpha_squared) * u_t_x + dt * u_dot;
    return u_tplus_x;
}

double wave(std::uint64_t t, std::uint64_t x)
{
    {
        std::lock_guard<hpx::mutex> l(u[t][x].mtx);
        //  hpx::util::format_to(cout, "calling wave... t={1} x={2}\n", t, x);
        if (u[t][x].computed)
        {
            //cout << ("already computed!\n");
            return u[t][x].u_value;
        }
        u[t][x].computed = true;

        if (t == 0)    //first timestep are initial values
        {
            //        hpx::util::format_to(cout, "first timestep\n");
            u[t][x].u_value = std::sin(
                2. * pi * static_cast<double>(x) * dx);    // initial u(x) value
            return u[t][x].u_value;
        }
    }

    // NOT using ghost zones here... just letting the stencil cross the periodic
    // boundary.
    future<double> n1;
    if (x == 0)
        n1 = async<wave_action>(here, t - 1, nx - 1);
    else
        n1 = async<wave_action>(here, t - 1, x - 1);

    future<double> n2 = async<wave_action>(here, t - 1, x);

    future<double> n3;
    if (x == (nx - 1))
        n3 = async<wave_action>(here, t - 1, 0);
    else
        n3 = async<wave_action>(here, t - 1, x + 1);

    double u_t_xminus = n1.get();    //get the futures
    double u_t_x = n2.get();
    double u_t_xplus = n3.get();

    if (t == 1)    //second time coordinate handled differently
    {
        std::lock_guard<hpx::mutex> l(u[t][x].mtx);
        double u_dot = 0;    // initial du/dt(x)
        u[t][x].u_value =
            calculate_u_tplus_x_1st(u_t_xplus, u_t_x, u_t_xminus, u_dot);
        return u[t][x].u_value;
    }
    else
    {
        double u_tminus_x = async<wave_action>(here, t - 2, x).get();
        std::lock_guard<hpx::mutex> l(u[t][x].mtx);
        u[t][x].u_value =
            calculate_u_tplus_x(u_t_xplus, u_t_x, u_t_xminus, u_tminus_x);
        return u[t][x].u_value;
    }
}

///////////////////////////////////////////////////////////////////////////////
int hpx_main(variables_map& vm)
{
    here = find_here();
    pi = 3.141592653589793238462643383279;

    //    dt = vm["dt-value"].as<double>();
    //    dx = vm["dx-value"].as<double>();
    //    c = vm["c-value"].as<double>();
    nx = vm["nx-value"].as<std::uint64_t>();
    nt = vm["nt-value"].as<std::uint64_t>();

    c = 1.0;

    dt = 1.0 / static_cast<double>(nt - 1);
    dx = 1.0 / static_cast<double>(nx - 1);
    alpha_squared = (c * dt / dx) * (c * dt / dx);

    // check that alpha_squared satisfies the stability condition
    if (0.25 < alpha_squared)
    {
        cout << (("alpha^2 = (c*dt/dx)^2 should be less than 0.25 for "
                  "stability!\n"))
             << std::flush;
    }

    u = std::vector<std::vector<data>>(nt, std::vector<data>(nx));

    hpx::util::format_to(cout, "dt = {1}\n", dt) << std::flush;
    hpx::util::format_to(cout, "dx = {1}\n", dx) << std::flush;
    hpx::util::format_to(cout, "alpha^2 = {1}\n", alpha_squared) << std::flush;

    {
        // Keep track of the time required to execute.
        high_resolution_timer t;

        std::vector<future<double>> futures;
        for (std::uint64_t i = 0; i < nx; i++)
            futures.push_back(async<wave_action>(here, nt - 1, i));

        // open file for output
        std::ofstream outfile;
        outfile.open("output.dat");

        auto f = [&](std::size_t i, hpx::future<double>&& f) {
            double x_here = static_cast<double>(i) * dx;
            hpx::util::format_to(outfile, "{1} {2}\n", x_here, f.get())
                << std::flush;
        };
        hpx::wait_each(f, futures);

        outfile.close();

        char const* fmt = "elapsed time: {1} [s]\n";
        hpx::util::format_to(std::cout, fmt, t.elapsed());
    }

    hpx::finalize();
    return 0;
}

///////////////////////////////////////////////////////////////////////////////
int main(int argc, char* argv[])
{
    // Configure application-specific options.
    options_description desc_commandline(
        "Usage: " HPX_APPLICATION_STRING " [options]");

    // clang-format off
    desc_commandline.add_options()
        ("dt-value",
           value<double>()->default_value(0.05),
           "dt parameter of the wave equation")

        ("dx-value",
           value<double>()->default_value(0.1),
           "dx parameter of the wave equation")

        ("c-value",
           value<double>()->default_value(1.0),
           "c parameter of the wave equation")

        ("nx-value",
           value<std::uint64_t>()->default_value(100),
           "nx parameter of the wave equation")

        ("nt-value",
           value<std::uint64_t>()->default_value(100),
           "nt parameter of the wave equation")
    ;
    // clang-format on

    // Initialize and run HPX.
    hpx::init_params init_args;
    init_args.desc_cmdline = desc_commandline;

    return hpx::init(argc, argv, init_args);
}
#endif

1	////////////////////////////////////////////////////////////////////////////////
2	// Copyright (c) 2012 Zach Byerly
3	// Copyright (c) 2011-2012 Bryce Adelstein-Lelbach
4	//
5	// SPDX-License-Identifier: BSL-1.0
6	// Distributed under the Boost Software License, Version 1.0. (See accompanying
7	// file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
8	////////////////////////////////////////////////////////////////////////////////
9
10	//
11	// This is a program written to evolve in time the equation:
12	//
13	// D^2 U / Dt^2 = c^2 D^2 U / Dx^2
14	//
15	// The parameter alpha = c*dt/dx must be less than 1 to ensure the stability
16	// of the algorithm.
17	// Discretizing the equation and solving for U(t+dt,x) yields
18	// alpha^2 * (U(t,x+dx)+U(t,x-dx))+2(1-alpha^2)*U(t,x) - U(t-dt,x)
19	//
20	// For the first timestep, we approximate U(t-dt,x) by u(t+dt,x)-2dtdu/dt(t,x)
21	//
22
23	// Include statements.
24	#include <hpx/config.hpp>
25	#if !defined(HPX_COMPUTE_DEVICE_CODE)
26	#include <hpx/chrono.hpp>
27	#include <hpx/hpx_init.hpp>
28	#include <hpx/include/async.hpp>
29	#include <hpx/iostream.hpp>
30	#include <hpx/modules/actions_base.hpp>
31	#include <hpx/modules/async_combinators.hpp>
32	#include <hpx/modules/format.hpp>
33
34	#include <cstddef>
35	#include <cstdint>
36	#include <fstream>
37	#include <iostream>
38	#include <mutex>
39	#include <vector>
40
41	using hpx::program_options::options_description;
42	using hpx::program_options::value;
43	using hpx::program_options::variables_map;
44
45	using hpx::id_type;
46	using hpx::invalid_id;
47
48	using hpx::async;
49	using hpx::future;
50
51	using hpx::chrono::high_resolution_timer;
52
53	using hpx::find_here;
54
55	using hpx::cout;
56
57	///////////////////////////////////////////////////////////////////////////////
58	// Globals.
59
60	//double const alpha_squared = 0.25;
61	double alpha_squared = 0;
62
63	// Initialized in hpx_main.
64	id_type here = invalid_id;
65	double pi = 0.;
66	double c = 0.;
67	double dt = 0.;
68	double dx = 0.;
69
70	// Command line argument.
71	std::uint64_t nt = 0;
72	std::uint64_t nx = 0;
73
74	struct data	10,100✔
75	{
76	// Default constructor: data d1;
77	data()
78	: mtx()	100✔
79	, u_value(0.0)	100✔
80	, computed(false)	100✔
81	{
82	}
83
84	// Copy constructor: data d1; data d2(d1);
85	// We can't copy the mutex, because mutexes are noncopyable.
86	data(data const& other)
87	: mtx()	10,000✔
88	, u_value(other.u_value)	10,000✔
89	, computed(other.computed)	10,000✔
90	{
91	}
92
93	data& operator=(data const& other)
94	{
95	u_value = other.u_value;
96	computed = other.computed;
97	return *this;
98	}
99
100	hpx::mutex mtx;
101	double u_value;
102	bool computed;
103	};
104
105	std::vector<std::vector<data>> u;
106
107	///////////////////////////////////////////////////////////////////////////////
108	// Forward declaration of the wave function.
109	double wave(std::uint64_t t, std::uint64_t x);
110
111	// Any global function needs to be wrapped into a plain_action if it should be
112	// invoked as a HPX-thread.
113	// This generates the required boilerplate we need for remote invocation.
114	HPX_PLAIN_ACTION(wave)	39,600✔
115
116	double calculate_u_tplus_x(	×
117	double u_t_xplus, double u_t_x, double u_t_xminus, double u_tminus_x)
118	{
119	double u_tplus_x = alpha_squared * (u_t_xplus + u_t_xminus) +	9,800✔
120	2.0 * (1 - alpha_squared) * u_t_x - u_tminus_x;	9,800✔
121	return u_tplus_x;	×
122	}
123
124	double calculate_u_tplus_x_1st(	×
125	double u_t_xplus, double u_t_x, double u_t_xminus, double u_dot)
126	{
127	double u_tplus_x = 0.5 * alpha_squared * (u_t_xplus + u_t_xminus) +	100✔
128	(1 - alpha_squared) * u_t_x + dt * u_dot;	100✔
129	return u_tplus_x;	×
130	}
131
132	double wave(std::uint64_t t, std::uint64_t x)	39,600✔
133	{
134	{
135	std::lock_guard<hpx::mutex> l(u[t][x].mtx);	39,600✔
136	// hpx::util::format_to(cout, "calling wave... t={1} x={2}\n", t, x);
137	if (u[t][x].computed)	39,600✔
138	{
139	//cout << ("already computed!\n");
140	return u[t][x].u_value;	29,600✔
141	}
142	u[t][x].computed = true;	10,000✔
143
144	if (t == 0) //first timestep are initial values	10,000✔
145	{
146	// hpx::util::format_to(cout, "first timestep\n");
147	u[t][x].u_value = std::sin(	100✔
148	2. * pi * static_cast<double>(x) * dx); // initial u(x) value	100✔
149	return u[t][x].u_value;	100✔
150	}
151	}
152
153	// NOT using ghost zones here... just letting the stencil cross the periodic
154	// boundary.
155	future<double> n1;	9,900✔
156	if (x == 0)	9,900✔
157	n1 = async<wave_action>(here, t - 1, nx - 1);	198✔
158	else
159	n1 = async<wave_action>(here, t - 1, x - 1);	19,602✔
160
161	future<double> n2 = async<wave_action>(here, t - 1, x);	9,900✔
162
163	future<double> n3;	9,900✔
164	if (x == (nx - 1))	9,900✔
165	n3 = async<wave_action>(here, t - 1, 0);	198✔
166	else
167	n3 = async<wave_action>(here, t - 1, x + 1);	19,602✔
168
169	double u_t_xminus = n1.get(); //get the futures	9,900✔
170	double u_t_x = n2.get();	9,900✔
171	double u_t_xplus = n3.get();	9,900✔
172
173	if (t == 1) //second time coordinate handled differently	9,900✔
174	{
175	std::lock_guard<hpx::mutex> l(u[t][x].mtx);	100✔
176	double u_dot = 0; // initial du/dt(x)
177	u[t][x].u_value =	100✔
178	calculate_u_tplus_x_1st(u_t_xplus, u_t_x, u_t_xminus, u_dot);
179	return u[t][x].u_value;
180	}
181	else
182	{
183	double u_tminus_x = async<wave_action>(here, t - 2, x).get();	29,400✔
184	std::lock_guard<hpx::mutex> l(u[t][x].mtx);	9,800✔
185	u[t][x].u_value =	9,800✔
186	calculate_u_tplus_x(u_t_xplus, u_t_x, u_t_xminus, u_tminus_x);
187	return u[t][x].u_value;
188	}
189	}
190
191	///////////////////////////////////////////////////////////////////////////////
192	int hpx_main(variables_map& vm)	1✔
193	{
194	here = find_here();	1✔
195	pi = 3.141592653589793238462643383279;	1✔
196
197	// dt = vm["dt-value"].as<double>();
198	// dx = vm["dx-value"].as<double>();
199	// c = vm["c-value"].as<double>();
200	nx = vm["nx-value"].as<std::uint64_t>();	2✔
201	nt = vm["nt-value"].as<std::uint64_t>();	2✔
202
203	c = 1.0;	1✔
204
205	dt = 1.0 / static_cast<double>(nt - 1);	1✔
206	dx = 1.0 / static_cast<double>(nx - 1);	1✔
207	alpha_squared = (c * dt / dx) * (c * dt / dx);	1✔
208
209	// check that alpha_squared satisfies the stability condition
210	if (0.25 < alpha_squared)	1✔
211	{
212	cout << (("alpha^2 = (c*dt/dx)^2 should be less than 0.25 for "
213	"stability!\n"))	1✔
214	<< std::flush;	1✔
215	}
216
217	u = std::vector<std::vector<data>>(nt, std::vector<data>(nx));	1✔
218
219	hpx::util::format_to(cout, "dt = {1}\n", dt) << std::flush;	1✔
220	hpx::util::format_to(cout, "dx = {1}\n", dx) << std::flush;	1✔
221	hpx::util::format_to(cout, "alpha^2 = {1}\n", alpha_squared) << std::flush;	2✔
222
223	{
224	// Keep track of the time required to execute.
225	high_resolution_timer t;
226
227	std::vector<future<double>> futures;	1✔
228	for (std::uint64_t i = 0; i < nx; i++)	101✔
229	futures.push_back(async<wave_action>(here, nt - 1, i));	200✔
230
231	// open file for output
232	std::ofstream outfile;	1✔
233	outfile.open("output.dat");	1✔
234
235	auto f = [&](std::size_t i, hpx::future<double>&& f) {	100✔
236	double x_here = static_cast<double>(i) * dx;	100✔
237	hpx::util::format_to(outfile, "{1} {2}\n", x_here, f.get())	100✔
238	<< std::flush;
239	};	1✔
240	hpx::wait_each(f, futures);
241
242	outfile.close();	1✔
243
244	char const* fmt = "elapsed time: {1} [s]\n";
245	hpx::util::format_to(std::cout, fmt, t.elapsed());	1✔
246	}	1✔
247
248	hpx::finalize();
249	return 0;	1✔
250	}
251
252	///////////////////////////////////////////////////////////////////////////////
253	int main(int argc, char* argv[])	1✔
254	{
255	// Configure application-specific options.
256	options_description desc_commandline(
257	"Usage: " HPX_APPLICATION_STRING " [options]");	1✔
258
259	// clang-format off
260	desc_commandline.add_options()	1✔
261	("dt-value",	1✔
262	value<double>()->default_value(0.05),	1✔
263	"dt parameter of the wave equation")
264
265	("dx-value",	1✔
266	value<double>()->default_value(0.1),	1✔
267	"dx parameter of the wave equation")
268
269	("c-value",	1✔
270	value<double>()->default_value(1.0),	1✔
271	"c parameter of the wave equation")
272
273	("nx-value",	1✔
274	value<std::uint64_t>()->default_value(100),	1✔
275	"nx parameter of the wave equation")
276
277	("nt-value",	1✔
278	value<std::uint64_t>()->default_value(100),	1✔
279	"nt parameter of the wave equation")
280	;
281	// clang-format on
282
283	// Initialize and run HPX.
284	hpx::init_params init_args;	1✔
285	init_args.desc_cmdline = desc_commandline;	1✔
286
287	return hpx::init(argc, argv, init_args);	1✔
288	}	1✔
289	#endif

STEllAR-GROUP / hpx / #882

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous