From 856e4e99bc87d99db01526232692c953cffee698 Mon Sep 17 00:00:00 2001 From: Michael Greenburg Date: Tue, 21 Nov 2023 13:50:43 -0700 Subject: [PATCH] Got rid of phase 1 example code that gave students a false sense of security --- project/overview.md | 95 ++++++--------------------------------------- project/phase1.md | 15 ++++--- project/phase9.md | 2 +- 3 files changed, 22 insertions(+), 90 deletions(-) diff --git a/project/overview.md b/project/overview.md index 222ec2b..afbb655 100644 --- a/project/overview.md +++ b/project/overview.md @@ -38,14 +38,14 @@ Here's how one step of `dt` for the whole wave rectangle might look in [Julia](. ```julia function step(u, v, dt) - m, n = size(u) + rows, cols = size(u) # Update v - for i in 2:m-1, j in 2:n-1 + for i in 2:rows-1, j in 2:cols-1 L = (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1]) / 2 - 2 * u[i,j] v[i,j] = (1 - dt * c) * v[i,j] + dt * L end # Update u - for i in 2:m-1, j in 2:n-1 + for i in 2:rows-1, j in 2:cols-1 u[i,j] += v[i,j] * dt end end @@ -69,17 +69,17 @@ Here's how the total energy might be determined in [Julia](../programming-resour ```julia function energy(u, v) - m, n = size(u) + rows, cols = size(u) E = 0 # Dynamic - for i in 2:m-1, j in 2:n-1 + for i in 2:rows-1, j in 2:cols-1 E += v[i,j]^2 / 2 end # Potential - for i in 1:m-1, j in 2:n-1 # along x axis (note i range) + for i in 1:rows-1, j in 2:cols-1 # along x axis (note i range) E += (u[i,j] - u[i+1,j])^2 / 4; end - for i in 2:m-1, j in 1:n-1 # along y axis (note j range) + for i in 2:rows-1, j in 1:cols-1 # along y axis (note j range) E += (u[i,j] - u[i,j+1])^2 / 4; end # Return the total @@ -91,7 +91,7 @@ end ## Running the Simulation -Your job is to determine an initial state (the specifics of which vary phase to phase), repeatedly step the simulation forward until total energy falls below an average of 0.001 per interior cell (i.e. **solve** the simulation), and do something with the resultant final state (again, the specifics vary by phase). Given damping coefficient `c`, time step `dt`, initial simulation time `t0`, initial displacement `u0`, and initial velocify `v0`, and functions [`step`](#moving-the-simulation-forward-in-time) and [`energy`](#stopping-criterion-energy) defined above, here is a Julia function that would solve the wave rectangle defined by `c`, `t0`, `u0`, and `v0` and return the changed elements of the final state. +Your job is to determine an initial state (the specifics of which vary phase to phase), repeatedly step the simulation forward until total energy falls below an average of 0.001 per interior cell (i.e. **solve** the simulation), and do something with the resultant final state (again, the specifics vary by phase). Given damping coefficient `c`, time step `dt`, initial simulation time `t0`, initial displacement `u0`, and initial velocity `v0`, and functions [`step`](#moving-the-simulation-forward-in-time) and [`energy`](#stopping-criterion-energy) defined above, here is a Julia function that would solve the wave rectangle defined by `c`, `t0`, `u0`, and `v0` and return the changed elements of the final state. ```julia function solve(c, dt, t0, u0, v0) @@ -99,8 +99,8 @@ function solve(c, dt, t0, u0, v0) t = t0 u, v = deepcopy.(u0, v0) # Constants - m, n = size(u) - stopping_energy = (m - 2) * (n - 2) / 1000 + rows, cols = size(u) + stopping_energy = (rows-2) * (cols-2) / 1000 # Solve while energy(u, v) > stopping_energy step(u, v, dt) @@ -123,80 +123,7 @@ Look at the example [code often](https://github.com/BYUHPC/sci-comp-course-examp -## Appendix B: Skeleton Wave Simulation Class - -Below is a class that contains the state of a wave simulation. It's provided since this isn't a programming class--we want you to learn relevant general concepts, not the painful ins and outs of C++. The parts you'll need to implement (for the [first phase](phase1.md), at least) are marked with `// TODO` comments. You'll likely need to modify it and split it into separate headers for each phase; look to the [example code](https://github.com/BYUHPC/sci-comp-course-example-cxx) for ideas on how to do so. - -```c++ -#include -#include -#include - -// Class representing a rectangular plane with Dirichlet boundary conditions over which waves propagate -// Called "orthotope" rather than "rectangle" since extra credit is given for generalizing to arbitrary dimension -class WaveOrthotope { - // Types - using size_type = size_t; - using value_type = double; - - // Members - const size_type m, n; // size - const value_type c; // damping coefficient - value_type t; // simulation time - std::vector disp, vel; // displacement and velocity arrays - -public: - // Constructor - // USAGE: auto my_wave_orthotope = WaveOrthotope(20, 30, 0.1, 0); - // SEE: https://www.learncpp.com/cpp-tutorial/constructor-member-initializer-lists/ - WaveOrthotope(auto m, auto n, auto c, value_type t=0): m(m), n(n), c(c), t(t), disp(m*n), vel(m*n) {} - - // Return the size as a std::array - // USAGE: const auto [rows, cols] = my_wave_orthotope.size(); - // SEE: https://codeburst.io/c-17-structural-binding-180696f7a678#63b1 - constexpr auto size() const { return std::array{m, n}; } - - // Displacement indexing - // USAGE: my_wave_orthotope.u(1, 2) = 1.0; - constexpr value_type u(auto i, auto j) const { return disp[i*n+j]; } - constexpr value_type& u(auto i, auto j) { return disp[i*n+j]; } - - // Displacement velocity indexing - // USAGE: const auto v12 = my_wave_orthotope.v(1, 2); - constexpr value_type v(auto i, auto j) const { return vel [i*n+j]; } - constexpr value_type& v(auto i, auto j) { return vel [i*n+j]; } - - // Return the energy contained in this WaveOrthotope - // USAGE: const auto E = my_wave_orthotope.energy(); - value_type energy() const { - value_type E{}; - // TODO: calculate total energy - return E; - } - - // Advance the membrane in time by dt - // USAGE: my_wave_orthotope.step(); - value_type step() { - const value_type dt = 0.01; - // TODO: update u and v by one time step - t += dt; - return t; - } - - // Advance the membrane in time by steps of dt until the average interior cell energy drops below 0.001 - // USAGE: const auto sim_time = my_wave_orthotope.solve(); - value_type solve() { - while (energy() > (m-2)*(n-2)*0.001) { - step(); - } - return t; - } -}; -``` - - - -## Appendix C: Mathematical Justification +## Appendix B: Mathematical Justification You don't need to read any of the following, but [some of it](#laplacian) is helpful for the extra credit. diff --git a/project/phase1.md b/project/phase1.md index c7bb790..c2ae465 100644 --- a/project/phase1.md +++ b/project/phase1.md @@ -39,17 +39,22 @@ g++ -std=c++20 -o wavesolve wavesolve.cpp You can test your implementation in the same way. You'll probably want to [set up VS Code](../resources.md#programming) on the supercomputer to develop your code unless you're familiar with vim or emacs. If you have no experience with [Linux](../lessons/2-linux.md), you can also [test with an online C++ compiler](../resources.md#compilation). -The [skeleton class in the overview](overview.md#appendix-b-skeleton-wave-simulation-class) and the `main` below constitute a good starting point, and I recommend using them as a scaffold for your program unless you have extensive experience with modern C++. There's not much need to worry about the [example code](https://github.com/BYUHPC/sci-comp-course-example-cxx) for this phase; it'll be more useful for structuring subsequent phases. +If you follow the [example code](https://github.com/BYUHPC/sci-comp-course-example-cxx) to create a `WaveOrthotope` class with `solve` and `sim_time` functions and a constructor taking rows, columns, and damping coefficient, your `main` can be very simple: ```c++ +#include +#include "whatever_your_header_is_called.hpp" + int main() { // Create WaveOrthotope - auto m = 25, n = 50; + auto rows = 25, cols = 50; auto c = 0.01; - auto w = WaveOrthotope(25, 50, 0.01); - // TODO: set interior cells of v to 0.1 + auto w = WaveOrthotope(rows, cols, c); + // Set interior cells of v to 0.1 + // TODO // Solve and print result - std::cout << w.solve() << std::endl; + w.solve() + std::cout << w.sim_time() << std::endl; return 0; } ``` diff --git a/project/phase9.md b/project/phase9.md index 6bb9f15..458882f 100644 --- a/project/phase9.md +++ b/project/phase9.md @@ -13,7 +13,7 @@ The self-contained [example code](https://github.com/BYUHPC/sci-comp-course-exam ### The Derivative Function -Defining the function that you'll pass to `SecondOrderODEProblem` itself isn't too bad, and `dhdt!` in the [example](https://github.com/BYUHPC/sci-comp-course-example-cxx/blob/main/src/initial.jl) offers some guidance. Notice that `dhdt!`'s main computation is just the [differential equation that defines the example problem](https://github.com/BYUHPC/sci-comp-course-example-cxx/tree/main#appendix-a-mathematical-justification)--you'll do something similar with the [damped wave equation differential equation](overview.md#appendix-c-mathematical-justification). Here's the gist: +Defining the function that you'll pass to `SecondOrderODEProblem` itself isn't too bad, and `dhdt!` in the [example](https://github.com/BYUHPC/sci-comp-course-example-cxx/blob/main/src/initial.jl) offers some guidance. Notice that `dhdt!`'s main computation is just the [differential equation that defines the example problem](https://github.com/BYUHPC/sci-comp-course-example-cxx/tree/main#appendix-a-mathematical-justification)--you'll do something similar with the [damped wave equation differential equation](overview.md#appendix-b-mathematical-justification). Here's the gist: ```julia function dvdt!(a, v, u, c, t) # u: displacement; v: velocity; a: acceleration