Code Base Review FY24

[psakievich]

Discussion to organize thoughts and comments while reviewing the code base. I'm logging things as I see them here so I can have a reference. I will clean up my thoughts and comments as I go.

Process:

1. Docs
2. Tests
3. Src

I want to see how well I understand the ideas from the tests first before looking at the actual implementations given the goal for TDD in the project plan.

@psakievich - Overview:

I see several positive things and potential for improvement. My high level summary is that the code's current form reads like a toolkit with a lot of fundamental math algorithms, and a few demonstration problems via unit-tests showing how you can piece the fundamental infrastructure together to solve a real problem.

It all seems like a good start. It is clear a lot of care and thought has gone into designing the low level data structures and algorithms.

I would like to see the code begin to take a general form through the main function for problem setup and things like BC/IC/Forcing terms to be applied generally. This would help to understand how the problems will be setup and used.
It will also help when considering what a coupling interface would look like for talking with other codes (openfast, nalu-wind, amr-wind etc).

Things I appreciate/like include:

• Variables are well named and it makes the code easier to read
• The low level math libraries are well tested and their use makes the algorithms easier to comprehend
• CI and CMake infrastructure are all looking good.
• Verification cases are built into the test suite and appear to be driving the development process
• Papers and references are in the code comments which make it easier to review and understand specific vernacular

Things I would like to see improved:

• There is a lot of allocations between methods that obscures the flow of the code (see my comment on void functions)
• Many of the tests are identical in function. Parameterization would help to comprehend the the extents of the testing better
• There is a lot of setup in the physics problem unit tests that would be nice to abstract into functions and methods.
• I see lots of hard coded values in the test suite. I would prefer to see them as variables unless there is a reference driving the value (i.e. from a paper). This will help with reduce duplicate code, as well as make the flow of the setup process and overall design.
• Documentation that shows the code working on device, convergence, verification results etc.

Things to anticipate and think about going forward

• Input deck file format and parser
• How abstract problems will be setup
  • For example your current physics problems in the test suite would be nice to refactor or reproduce with a minimum input deck that you either hard code in a header or add as a first round of regression tests.
• API's for the abstract problems via input files and code to code communication (coupling with Nalu-Wind for example)
• Restarts, and interactions with a coupled aero solver via non-linear iterations
• Do we want to support variable time stepping? Can things be dynamically changed upon restart?

[psakievich]

openturbine/docs/compiling.md

Lines 13 to 15 in 45aab4d

	when available. For workstations, package managers such as
	conda, APT, and Homebrew will provide the easiest experience
	for linking and runtime path searches.

Probably worth getting started with spack on this project since the rest of the exawind suite is based around using spack.

[psakievich]

openturbine/tests/unit_tests/gen_alpha_poc/test_quaternions.cpp

Lines 219 to 234 in 45aab4d

	ASSERT_NEAR(v.GetXComponent(), expected.GetXComponent(), 1e-20);
	ASSERT_NEAR(v.GetYComponent(), expected.GetYComponent(), 1e-20);
	ASSERT_NEAR(v.GetZComponent(), expected.GetZComponent(), 1e-20);
	}
	TEST(QuaternionTest, QuaternionFromAngleAxis_ZeroAngle) {
	double angle = 0.;
	Vector axis{1., 0., 0.};
	Quaternion q = quaternion_from_angle_axis(angle, axis);
	Quaternion expected(1., 0., 0., 0.);
	ASSERT_NEAR(q.GetScalarComponent(), expected.GetScalarComponent(), 1e-6);
	ASSERT_NEAR(q.GetXComponent(), expected.GetXComponent(), 1e-6);
	ASSERT_NEAR(q.GetYComponent(), expected.GetYComponent(), 1e-6);
	ASSERT_NEAR(q.GetZComponent(), expected.GetZComponent(), 1e-6);

Pretty wide range of tolerances here. I am curious as to why, and if there are separate tolerances expected for single/double precision mode. I would suggest making all similar tolerances a variable, and using something like constexpr if to set it appropriately based on precision/architecture.

I suspect that really understanding and monitoring the tolerances for these basic math operations will be important for everything else going forward.

[psakievich]

Similar question in this file. Are you building in the flexibility for multiple platforms with ranging tolerances etc? Point of discussion, perhaps this precision is sufficient to satisfy the testing already.

openturbine/tests/unit_tests/gen_alpha_poc/test_heavy_top.cpp

Lines 92 to 102 in 45aab4d

	expect_kokkos_view_1D_equal(
	residual_vector,
	{
	14., // row 1
	13., // row 2
	159.150000, // row 3
	3.135087, // row 4
	0.46875, // row 5
	-1.310708 // row 6
	}
	);

[faisal-bhuiyan]

Recorded in #123.

[alanw0]

Low-level comment, initialization of views is sometimes a significant performance issue.

openturbine/src/gen_alpha_poc/generalized_alpha_time_integrator.cpp

Lines 143 to 147 in 45aab4d

	// Precondition the linear solve (Bottasso et al 2008)
	const auto dl = Kokkos::View<double**>("dl", size + n_constraints, size + n_constraints);
	const auto dr = Kokkos::View<double**>("dr", size + n_constraints, size + n_constraints);
	Kokkos::deep_copy(dl, 0.);
	Kokkos::deep_copy(dr, 0.);

I think that Kokkos Views are zero-initialized by default, possibly meaning that the deep_copy with zero is un-necessary.
In fact a common performance optimization (sometimes significant) is to make sure Kokkos::View doesn't initialize if you know that you will be writing into it (not +=) and don't care what the initial value is, such as here:

openturbine/src/gen_alpha_poc/generalized_alpha_time_integrator.cpp

Lines 201 to 202 in 45aab4d

	auto soln_increments = Kokkos::View<double*>("soln_increments", residuals.size());
	Kokkos::deep_copy(soln_increments, residuals);

In this case you would declare the Kokkos::View with the Kokkos::WithoutInitializing argument to avoid the internal kernel that is launched to perform the initialization.

[faisal-bhuiyan]

Thanks for the suggestion, @alanw0. Recorded as an issue in the repo: #118.

[alanw0]

Regarding passing 'raw' Kokkos::View types everywhere:

openturbine/src/gen_alpha_poc/linearization_parameters.h

Lines 14 to 18 in 45aab4d

	virtual Kokkos::View<double*> ResidualVector(
	const Kokkos::View<double*> /* gen_coords */, const Kokkos::View<double*> /* velocity */,
	const Kokkos::View<double*> /* acceleration */,
	const Kokkos::View<double*> /* lagrange_multipliers */
	) = 0;

It's nice to see what the actual type is, but if performance optimizations or development leads to adding extra template parameters such as textured-memory (Kokkos::MemoryRandomAccess), or Atomic traits, then we may wish for central typedefs so that changes can be minimized. (These kinds of changes would manifest as extra template parameters in the Kokkos::View type.)

[deslaughter]

I think named types would be helpful since many of these have been refactored into multidimensional arrays where some of the dimensions are fixed size. It would remove the need to specify the dimension size in the type wherever it's being used.

[alanw0]

minor comment: should arguments be named in header declarations?

openturbine/src/gen_alpha_poc/time_integrator.h

Lines 21 to 22 in 45aab4d

	virtual std::vector<State>
	Integrate(const State&, size_t, std::shared_ptr<LinearizationParameters>) = 0;

I don't know what the second argument is (size_t) unless I go find an override of this method in a derived class. It turns out (e.g. in generalized_alpha_time_integrator.cpp) that the size_t argument is n_constraints.

In general: the compiler doesn't require argument names for declarations, only types. I like argument names, although I recognize there is a "don't repeat yourself" viewpoint which says to only name the arguments in the cpp file where it is needed. I would suggest to just be consistent with the chosen style. There are 3 possible styles:

void foo(int);
void foo(int /*argname*/);
void foo(int argname);

I've seen all 3 of these styles in the Open Turbine code base.

[faisal-bhuiyan]

I don't know what the second argument is (size_t) unless I go find an override of this method in a derived class. It turns out (e.g. in generalized_alpha_time_integrator.cpp) that the size_t argument is n_constraints.

I agree with this observation - I find argument names to be helpful for primitive types.

I would suggest to just be consistent with the chosen style.

Looks like we need to nail down the style choice in OpenTurbine - will start a style guide discussion to handle this.

[faisal-bhuiyan]

Transferred as #122

[alanw0]

The large number of unit tests is great, this is a nice artifact of doing test driven development. As development proceeds towards having an API for the "entry point" where a calling application would use Open Turbine, I would recommend creating more "doc" tests in addition to unit tests. Unit tests often check many aspects of correctness and "corner cases", and may not be perfectly readable for all audiences.
Doc tests should not add additional test coverage but rather act as a demonstration or example of how a client should interact with the library.
This test seems close to what I might consider a doc test:

openturbine/tests/unit_tests/gen_alpha_poc/test_generalized_alpha_solver.cpp

Line 58 in 45aab4d

TEST(TimeIntegratorTest, AdvanceAnalysisTimeByNumberOfSteps) {

[faisal-bhuiyan]

@alanw0 Copy the above, will try to keep the (additional) utility of unit tests to work as a documentation for the codebase in mind. Please feel free to share some thoughts on the details and an example or two (in a codebase), if possible.

[psakievich]

Minor nitpick: some of these tests would be good candidates for a parameterization refactor.
Example:

openturbine/tests/unit_tests/gen_alpha_poc/test_math_utilities.cpp

Lines 8 to 30 in 45aab4d

	TEST(MathUtilitiesTest, CloseTo_Set1) {
	ASSERT_TRUE(close_to(1., 1.));
	}
	TEST(MathUtilitiesTest, CloseTo_Set2) {
	ASSERT_TRUE(close_to(1., 1. + kTolerance / 10.));
	}
	TEST(MathUtilitiesTest, CloseTo_Set3) {
	ASSERT_TRUE(close_to(1., 1. - kTolerance / 10.));
	}
	TEST(MathUtilitiesTest, CloseTo_Set4) {
	ASSERT_FALSE(close_to(1., 1. + kTolerance * 10.));
	}
	TEST(MathUtilitiesTest, CloseTo_Set5) {
	ASSERT_FALSE(close_to(1., 1. - kTolerance * 10.));
	}
	TEST(MathUtilitiesTest, CloseTo_Set6) {
	ASSERT_TRUE(close_to(kTolerance / 10., kTolerance / 10.));
	}

[faisal-bhuiyan]

Totally agree with the above - will create an issue to take care of this.

[faisal-bhuiyan]

Recorded here in #115.

[psakievich]

GLL quadrature is used a lot in the tests. I would just make a class for it. Seems like that is what you will use in the actual end application once you set it up as well.

openturbine/src/gebt_poc/quadrature.h

Line 23 in 45aab4d

class UserDefinedQuadrature : public Quadrature {

[michaelasprague]

As a point of clarification, we will use other quadrature methods as well (most likely trapezoidal when we have variable properties).

[faisal-bhuiyan]

I agree - we should add the GLL quadrature as a derived class from Quadrature and use that in the test suite to avoid repetition.

[faisal-bhuiyan]

Follow up issue: #116

[psakievich]

openturbine/src/gebt_poc/solver.cpp

Lines 172 to 264 in 45aab4d

	void ElementalStaticForcesResidual(
	const Kokkos::View<double*> position_vectors, const Kokkos::View<double*> gen_coords,
	const StiffnessMatrix& stiffness, const Quadrature& quadrature, Kokkos::View<double*> residual
	) {
	const auto n_nodes = gen_coords.extent(0) / kNumberOfLieAlgebraComponents;
	const auto order = n_nodes - 1;
	const auto n_quad_pts = quadrature.GetNumberOfQuadraturePoints();
	auto nodes = Kokkos::View<double* [3]>("nodes", n_nodes);
	for (std::size_t i = 0; i < n_nodes; ++i) {
	auto index = i * kNumberOfLieAlgebraComponents;
	Kokkos::deep_copy(
	Kokkos::subview(nodes, i, Kokkos::ALL),
	Kokkos::subview(position_vectors, Kokkos::make_pair(index, index + 3))
	);
	}
	// Allocate Views for some required intermediate variables
	auto gen_coords_qp = Kokkos::View<double[kNumberOfLieAlgebraComponents]>("gen_coords_qp");
	auto gen_coords_derivatives_qp =
	Kokkos::View<double[kNumberOfLieAlgebraComponents]>("gen_coords_derivatives_qp");
	auto position_vector_qp =
	Kokkos::View<double[kNumberOfLieAlgebraComponents]>("position_vector_qp");
	auto pos_vector_derivatives_qp =
	Kokkos::View<double[kNumberOfLieAlgebraComponents]>("pos_vector_derivatives_qp");
	auto curvature = Kokkos::View<double[kNumberOfVectorComponents]>("curvature");
	auto sectional_strain = Kokkos::View<double[kNumberOfLieGroupComponents]>("sectional_strain");
	auto sectional_stiffness =
	Kokkos::View<double[kNumberOfLieGroupComponents][kNumberOfLieGroupComponents]>(
	"sectional_stiffness"
	);
	Kokkos::deep_copy(residual, 0.);
	for (size_t i = 0; i < n_nodes; ++i) {
	const auto node_count = i;
	for (size_t j = 0; j < n_quad_pts; ++j) {
	// Calculate required interpolated values at the quadrature point
	const auto q_pt = quadrature.GetQuadraturePoints()[j];
	auto shape_function = LagrangePolynomial(order, q_pt);
	auto shape_function_derivative = LagrangePolynomialDerivative(order, q_pt);
	auto shape_function_vector = gen_alpha_solver::create_vector(shape_function);
	auto shape_function_derivative_vector =
	gen_alpha_solver::create_vector(shape_function_derivative);
	auto jacobian = CalculateJacobian(nodes, shape_function_derivative_vector);
	InterpolateNodalValues(gen_coords, shape_function, gen_coords_qp);
	InterpolateNodalValueDerivatives(
	gen_coords, shape_function_derivative, jacobian, gen_coords_derivatives_qp
	);
	InterpolateNodalValues(position_vectors, shape_function, position_vector_qp);
	InterpolateNodalValueDerivatives(
	position_vectors, shape_function_derivative, jacobian, pos_vector_derivatives_qp
	);
	// Calculate the curvature and sectional strain
	NodalCurvature(gen_coords_qp, gen_coords_derivatives_qp, curvature);
	CalculateSectionalStrain(
	pos_vector_derivatives_qp, gen_coords_derivatives_qp, curvature, sectional_strain
	);
	// Calculate the sectional stiffness matrix in inertial basis
	auto rotation_0 = gen_alpha_solver::EulerParameterToRotationMatrix(
	Kokkos::subview(position_vector_qp, Kokkos::make_pair(3, 7))
	);
	auto rotation = gen_alpha_solver::EulerParameterToRotationMatrix(
	Kokkos::subview(gen_coords_qp, Kokkos::make_pair(3, 7))
	);
	SectionalStiffness(stiffness, rotation_0, rotation, sectional_stiffness);
	// Calculate elastic forces i.e. F^C and F^D vectors
	auto elastic_forces_fc =
	Kokkos::View<double[kNumberOfLieGroupComponents]>("elastic_forces_fc");
	auto elastic_forces_fd =
	Kokkos::View<double[kNumberOfLieGroupComponents]>("elastic_forces_fd");
	NodalElasticForces(
	sectional_strain, rotation, pos_vector_derivatives_qp, gen_coords_derivatives_qp,
	sectional_stiffness, elastic_forces_fc, elastic_forces_fd
	);
	// Calculate the residual at the quadrature point
	const auto q_weight = quadrature.GetQuadratureWeights()[j];
	Kokkos::parallel_for(
	kNumberOfLieGroupComponents,
	KOKKOS_LAMBDA(const size_t i) {
	residual(node_count * kNumberOfLieGroupComponents + i) +=
	q_weight *
	(shape_function_derivative_vector(node_count) * elastic_forces_fc(i) +
	jacobian * shape_function_vector(node_count) * elastic_forces_fd(i));
	}
	);
	}
	}
	}

If you can get away with not using void functions to update parameters it really helps with readability. This pattern is repeated all throughout the code. It is hard to know what is getting updated and what isn't in these types of functions.

In the instance above, I would prefer a return of the residual in this case so it is clear what the output is. Also, in this case you are zero'ing out the residual anyways so it doesn't seem like you need to pass by reference. There is an argument for saving on memory allocations, but looking at the utilization of these and the memory is getting allocated right before calling it:

openturbine/src/gebt_poc/clamped_beam.cpp

Lines 75 to 78 in 45aab4d

	auto residual_gen_coords = Kokkos::subview(residual, Kokkos::make_pair(zero, size_dofs));
	ElementalStaticForcesResidual(
	position_vectors_, gen_coords_1D, stiffness_matrix_, quadrature_, residual_gen_coords
	);

Which I would prefer to see as:

    auto residual_gen_coords = ElementalStaticForcesResidual(
        position_vectors_, gen_coords_1D, stiffness_matrix_, quadrature_);

[faisal-bhuiyan]

I agree with the above, returning a type from a function is helpful as a reader to understand what's being done in the function. I can take care of it as part of the refactor.

@deslaughter Could you remind me of your point regarding the above? I remember we discussed the above issue in our 2/20 meeting but cannot recall what it was.

[faisal-bhuiyan]

Referenced in the following issue: #117.

[psakievich]

Group Discussion: @rafmudaf @deslaughter @psakievich @alanw0

Languages/Core Libraries (C++/Kokkos/Julia):

• Overall none of the reviewers see a compelling reason to go to Julia. The cost of adding in an entirely new ecosystems of build dependencies seems like a high cost
• Kokkos Kernels were difficult for debugging and lacked a lot of linear algebra features @deslaughter
  • Prototyped in rust and that seems viable. There are GPU libraries available, but we do not have active expertise there.
• Rather than looking for a new technology perhaps we should focus on using the libraries and technology we already have expertise in well/better @rafmudaf. Specifically looking to Kokkos+Hypre(for sparse linear systems)
• Linear algebra: seems like what we need is a clean API to translate the matrices to make it easier to port in and test various solver libraries.

[deslaughter]

I would like to keep the Rust development gonig as I've found it immensely helping for debugging the theory and prototyping the code. I also think it could serve as a good comparison for performance, in preliminary testing it's already slightly faster than the BeamDyn implementation in OpenFAST. As mentioned above, there are linear algebra libraries available such as ArrayFire which includes GPU acceleration and sparse matrix operations (also has C++ and Python APIs).

[psakievich]

This seems reasonable to me. The main concern I would have is the added infrastructure and context switching a new language and libraries provide. So I would just say keep that in mind. If the rust implementation proves to be significantly easier then I would ask to work out how to deploy it via spack on multiple machines before committing to a switch down the road.

[alanw0]

Question from Mike: Thoughts on the overall structure/organization of the code.

[psakievich]

The organization seems mostly logical, but its current organization seems more like a tool kit than an application. Specifically, the unit-tests show how concrete problems can be solved by assembling the tool kit components, but the abstraction of a complete problem (i.e read data, create objects, apply bc's forces etc, solve) is missing.

I would like to see the code reach a level of maturity where an arbitrary problem can be specified with a minimal set of information. Something like:

problem:
 name: "clamped beam"
 beam_positions:
 - start_position: 0, 0, 0
   end_position: 1, 0, 0
 boundary_conditions:
 - blah
 time_integration:
 - start_time: 0
 - end_time: 1
 - time_step: 1e-3
 solver_settings:
   type: kokkos_kernels # could be expanded to something else
   tolerances:
      absolute: 1e-12
      relative: 1e-7
   preconditioner_params:
   - doe

The details above are not important, or rather not meant to be a precise description, more that the idea of what is the minimal working set of data to describe the problem has been explored and captured in a clean form. This will allow for things like solver abstractions (adding hypre, trilinos, petsc etc), different quadrature types.

Once a first draft of this is in place I think we will have a lot more to discuss regarding structure.

I would encourage the team to write out examples of what they would want to put into an input file to describe these problems and then explore what it would take to turn that into code.
Anything that can be a parameter should be, but also utilize defaults so the final form of the input deck (and eventual API) is not overwhelming and pedantic.

Realistically parameterizing everything is a lot of work though. If you design your parsing and data structures (use structs etc) to allow for backward compatible additions then you don't need to define all the parameters at once. You can mark them as parametrizable in the code with a comment or something and turn them into a parameter with the default being set to what the hard coded option was to preserve backwards compatibility. Omission of new parameters in the input deck and/or API preserves old behavior.

[deslaughter]

A significant refactoring of the code is needed to minimize the dynamic allocation of Kokkos Views. Since the size of the problem doesn't change after initialization, all of the system matrices could be allocated once and subviews passed to relevant functions. The GEBT interpolation (shape function) arrays can be calculated at initialization and stored in matrices which can then interpolate the coordinates with matrix multiplication instead of individual operations. The Rust prototype was structured to minimize allocations and may be a good reference for the matrix sizing and operations. I think some performance profiling would be helpful to see which parts are being called repeatedly that may not need to be.

I agree with @psakievich's comment that the code seems more like a toolkit than a program. Do we want to have the ability to use it as a toolkit with the main program being geared towards modeling a turbine? This approach may help us keep the elements and solver more general and extensible. Even if we don't, one of the next steps will need to be defining the representation of the turbine model in OpenTurbine, which could then guide how the element and solver API is structured.

Personally, I'd like to clean up the existing code and generalize the element and solver APIs so that the turbine model can be constructed in a generic manor. It has been difficult to represent more complex turbine designs (multi-rotor) in OpenFAST because of its design assumptions.

[michaelasprague]

@deslaughter is the refactor to minimize dynamic allocation something you are interested in undertaking?

Regarding the toolkit aspects, I like that the work so far is viewed as a toolkit -- we can now add the layer(s) necessary to bring the parts together for a HAWT simulation -- but I don't want hinder our ability to build other systems.

[deslaughter]

Yes, I'd like to take a stab at doing the refactorization. I'm not sure when I can fit it in, let's discuss next week.

[alanw0]

Question from Mike: Thoughts on the code's readiness for an API for coupling to other codes.

[psakievich]

I detailed a lot of thoughts on this here, but I will expand on the idea related to coupling.

The ability to express a minimal problem is critical for a successful API. The steps I see are:

1. Identify the minimum input data needed for openturbine to do it's work
2. Assemble data structures that can be included via headers to express this information in simple types
3. Make this headers into a header only library that both codes can include

Until the data is cleanly assembled for running arbitrary classes of problems on the openturbine side, it is probably too early to be discuss much more on the API.

[rafmudaf]

I considered OpenTurbine from a high level perspective. After our group discussion (summarized in #108 (comment)), I got more insight into the current state of OpenTurbine.

In particular, the architecture of the software as a wind turbine analysis tool has not yet been developed. My primary feedback was similar to #108 (reply in thread) in that there are currently the building blocks for a cohesive software but no clear architecture yet. It sounds like now is the time to begin addressing this, and I recommend to step back and consider the design of this from a top-down, holistic perspective.

I also noticed there is a particular coding style (i.e. the trailing _ on properties, the GetXYZ methods, the place where long lines wrap, etc), but I don't see this style documented. As the number of developers scale, this will become more important. I suggest to start a style guide discussion to capture these low-level decisions with the understanding that it will change over time. See NREL/floris#292 for an idea of how this could be structured.

Is doxygen configured? I see the doxygen directory, but I don't see an entry point in the documentation site.

Similar to #108 (comment), I suggest to add some reporting on the problems you've solved so far. I would treat this as verification documentation, so add a static page to the docs describing the problem, how OpenTurbine was extended to solve the problem, and show the solution compared to any analytical results. Include a mesh dependency and convergence analysis.

[faisal-bhuiyan]

I also noticed there is a particular coding style (i.e. the trailing _ on properties, the GetXYZ methods, the place where long lines wrap, etc), but I don't see this style documented. As the number of developers scale, this will become more important. I suggest to start a style guide discussion to capture these low-level decisions with the understanding that it will change over time. See NREL/floris#292 for an idea of how this could be structured.

Thanks for the suggestion - I agree the style choices might appear random to the reader. Started a discussion on this here:
#125