From 01cb8e1b01190c55fac7c6fdfa3d44811a55ee0a Mon Sep 17 00:00:00 2001 From: Matt Rolchigo Date: Fri, 12 Sep 2025 16:28:56 -0400 Subject: [PATCH] Separated steering vector construction for active, other cells --- bin/run.cpp | 2 +- bin/runCoupled.cpp | 4 +- src/CAcelldata.hpp | 107 ++--- src/CAnucleation.hpp | 108 ++++- src/CAprint.hpp | 40 +- src/CAtemperature.hpp | 796 +++++++++++++++++------------------ src/CAtimers.hpp | 12 +- src/CAupdate.hpp | 389 +++++++++-------- src/runCA.hpp | 60 +-- unit_test/tstCellData.hpp | 76 +--- unit_test/tstInterface.hpp | 187 ++++---- unit_test/tstNucleation.hpp | 111 ++--- unit_test/tstTemperature.hpp | 86 ++-- unit_test/tstUpdate.hpp | 7 +- 14 files changed, 1032 insertions(+), 953 deletions(-) diff --git a/bin/run.cpp b/bin/run.cpp index fa7c4f64..5cdb1463 100644 --- a/bin/run.cpp +++ b/bin/run.cpp @@ -50,7 +50,7 @@ int main(int argc, char *argv[]) { Grid grid(inputs.simulation_type, id, np, inputs.domain.number_of_layers, inputs.domain, inputs.substrate, inputs.temperature); // Temperature fields characterized by data in this structure - Temperature temperature(grid, inputs.temperature, inputs.print.store_solidification_start); + Temperature temperature(grid, inputs.temperature, inputs.print); runExaCA(id, np, inputs, timers, grid, temperature); } diff --git a/bin/runCoupled.cpp b/bin/runCoupled.cpp index 49c1cc3d..21ccbba2 100644 --- a/bin/runCoupled.cpp +++ b/bin/runCoupled.cpp @@ -60,8 +60,8 @@ void runCoupled(int id, int np, Finch::Inputs finch_inputs, Inputs exaca_inputs) Grid exaca_grid(id, np, exaca_inputs.domain, exaca_inputs.substrate, exaca_inputs.temperature, finch_inputs.space.cell_size, exaca_low_corner, exaca_high_corner); // Temperature fields characterized by data in this structure - Temperature temperature(id, np, exaca_grid, exaca_inputs.temperature, app.getSolidificationData(), - exaca_inputs.print.store_solidification_start); + Temperature temperature(id, np, exaca_grid, exaca_inputs.temperature, exaca_inputs.print, + app.getSolidificationData()); // Now run ExaCA runExaCA(id, np, exaca_inputs, timers, exaca_grid, temperature); diff --git a/src/CAcelldata.hpp b/src/CAcelldata.hpp index d6a18c0a..ce78a4d2 100644 --- a/src/CAcelldata.hpp +++ b/src/CAcelldata.hpp @@ -10,6 +10,7 @@ #include "CAgrid.hpp" #include "CAinputdata.hpp" #include "CAparsefiles.hpp" +#include "CAtemperature.hpp" #include "CAtypes.hpp" #include "mpi.h" @@ -41,10 +42,8 @@ struct CellData { int next_layer_first_epitaxial_grain_id = 1; int num_prior_nuclei = 0; - view_type_int grain_id_all_layers; + view_type_int grain_id_all_layers, cell_type; view_type_short phase_id_all_layers; - view_type_int cell_type; - view_type_short layer_id_all_layers; view_type_bool melt_edge, melt_edge_all_layers; // Substrate inputs from file SubstrateInputs _inputs; @@ -56,10 +55,8 @@ struct CellData { // cell_type only exists for the current layer of a multilayer problem CellData(const Grid &grid, SubstrateInputs inputs, const bool store_melt_pool_edge = false) : grain_id_all_layers(view_type_int("GrainID", grid.domain_size_all_layers)) - , phase_id_all_layers(view_type_short("PhaseID", grid.domain_size_all_layers)) , cell_type(view_type_int(Kokkos::ViewAllocateWithoutInitializing("cell_type"), grid.domain_size)) - , layer_id_all_layers( - view_type_short(Kokkos::ViewAllocateWithoutInitializing("LayerID"), grid.domain_size_all_layers)) + , phase_id_all_layers(view_type_short("PhaseID", grid.domain_size_all_layers)) , _inputs(inputs) , _store_melt_pool_edge(store_melt_pool_edge) { if (_store_melt_pool_edge) { @@ -77,7 +74,7 @@ struct CellData { // Initializes the single active cell and associated active cell data structures for the single grain at the domain // center - void initSubstrate(const int id, const Grid &grid) { + void initSubstrate_SingleGrain(const int id, const Grid &grid) { // Location of the single grain int grain_location_x = Kokkos::floorf(static_cast(grid.nx) / 2.0); @@ -98,7 +95,7 @@ struct CellData { int coord_y_global = coord_y + grid.y_offset; if ((coord_x == grain_location_x) && (coord_y_global == grain_location_y) && (coord_z == grain_location_z)) { - cell_type_local(index) = FutureActive; + cell_type_local(index) = Active; grain_id_all_layers_local(index) = single_grain_orientation_local + 1; } else @@ -189,7 +186,7 @@ struct CellData { // Initializes cell types and epitaxial Grain ID values where substrate grains are future active cells on the bottom // surface of the constrained domain - void initSubstrate(const int id, const Grid &grid, const unsigned long rng_seed) { + void initSubstrate_Directional(const int id, const Grid &grid, const unsigned long rng_seed) { // Fill the view of cell X, Y, and ID values, updating the number of substrate active cells appropriately // TODO: Could generate random numbers on GPU, instead of using host view and copying over - but would also need @@ -200,9 +197,7 @@ struct CellData { auto act_cell_data = Kokkos::create_mirror_view_and_copy(memory_space(), act_cell_data_host); // Start with all cells as liquid prior to locating substrate grain seeds - // All cells have LayerID = 0 as this is not a multilayer problem Kokkos::deep_copy(cell_type, Liquid); - Kokkos::deep_copy(layer_id_all_layers, 0); // Local copies for lambda capture. auto cell_type_local = cell_type; @@ -221,7 +216,7 @@ struct CellData { int coord_y = act_cell_data(n, 1) - grid.y_offset; int coord_z = 0; int index = grid.get1DIndex(coord_x, coord_y, coord_z); - cell_type_local(index) = FutureActive; + cell_type_local(index) = Active; grain_id_all_layers_local(index) = act_cell_data(n, 2); // assign GrainID > 0 to epitaxial seeds } }); @@ -262,7 +257,7 @@ struct CellData { // These cells are also future active cells // For directional solidification, only one "layer" so index for all layers can be used here to // index cell_type - cell_type_local(index_all_layers) = FutureActive; + cell_type_local(index_all_layers) = Active; } }); } @@ -275,8 +270,7 @@ struct CellData { // Initializes substrate grain structure using either a Voronoi assignment of grain ID values to cells or reading // the data from a file - void initSubstrate(const int id, const Grid &grid, const unsigned long rng_seed, - view_type_int number_of_solidification_events) { + void initSubstrate_BaseplatePowder(const int id, const Grid &grid, const unsigned long rng_seed) { // Determine the number of cells in the Z direction that are part of the baseplate int baseplate_size_z = getBaseplateSizeZ(id, grid); @@ -309,12 +303,9 @@ struct CellData { std::cout << "Expected number of grains in the powder layers of this simulation: " << number_total_powder_grains << std::endl; - // LayerID starts at -1 for all cells - Kokkos::deep_copy(layer_id_all_layers, -1); - - // Initialize cell types and layer IDs based on whether cells will solidify in layer 0 or not - initCellTypeLayerID(0, id, grid, number_of_solidification_events); - + // Resize view and init cell types to solid + Kokkos::realloc(cell_type, grid.domain_size); + Kokkos::deep_copy(cell_type, Solid); MPI_Barrier(MPI_COMM_WORLD); if (id == 0) { std::cout << "Grain struct initialized" << std::endl; @@ -633,8 +624,7 @@ struct CellData { } // Sets up views, powder layer (if necessary), and cell types for the next layer of a multilayer problem - void initNextLayer(const int nextlayernumber, const int id, const Grid &grid, const unsigned long rng_seed, - view_type_int number_of_solidification_events) { + void initNextLayer(const int nextlayernumber, const int id, const Grid &grid, const unsigned long rng_seed) { // Subviews for the next layer's grain id, layer id, cell type are constructed based on updated layer bound // z_layer_bottom @@ -650,50 +640,17 @@ struct CellData { powder_rng_seed, powder_layer_label); } - // Initialize active cell data structures and nuclei locations for the next layer "layernumber + 1" - initCellTypeLayerID(nextlayernumber, id, grid, number_of_solidification_events); + // Resize and init cell types to Solid for layer "layernumber + 1" + Kokkos::realloc(cell_type, grid.domain_size); + Kokkos::deep_copy(cell_type, Solid); // Current layer melt pool edge indicator, if needed if (_store_melt_pool_edge) getCurrentLayerMeltEdge(grid.layer_range); } - // Initializes cells for the current layer as either solid (don't resolidify) or tempsolid (will melt and - // resolidify) - void initCellTypeLayerID(const int layernumber, const int id, const Grid &grid, - view_type_int number_of_solidification_events) { - - int melt_pool_cell_count; - // Realloc celltype to the domain size of the next layer - Kokkos::realloc(cell_type, grid.domain_size); - // Local copies for lambda capture. - auto cell_type_local = cell_type; - auto layer_id_all_layers_local = layer_id_all_layers; - - auto policy = Kokkos::RangePolicy(0, grid.domain_size); - - Kokkos::parallel_reduce( - "cell_typeInitSolidRM", policy, - KOKKOS_LAMBDA(const int &index, int &local_count) { - int index_all_layers = index + grid.z_layer_bottom * grid.nx * grid.ny_local; - if (number_of_solidification_events(index) > 0) { - cell_type_local(index) = TempSolid; - layer_id_all_layers_local(index_all_layers) = layernumber; - local_count++; - } - else - cell_type_local(index) = Solid; - }, - melt_pool_cell_count); - int total_melt_pool_cell_count; - MPI_Reduce(&melt_pool_cell_count, &total_melt_pool_cell_count, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); - if (id == 0) - std::cout << "Number of cells across all ranks to undergo solidification at least once: " - << total_melt_pool_cell_count << std::endl; - } - // Stores/returns the volume fraction of nucleated grains to the console // Moved from CAfunctions.hpp - float calcVolFractionNucleated(const int id, const Grid &grid) { + float calcVolFractionNucleated(const int id, const Grid &grid, Temperature &temperature) { // For interior cells, add the number of cells that underwent melting/solidification and the number of cells // with sub-zero grain IDs @@ -701,7 +658,6 @@ struct CellData { int nucleated_grain_cells_local = 0; // Local copies for lambda capture. auto grain_id_all_layers_local = grain_id_all_layers; - auto layer_id_all_layers_local = layer_id_all_layers; Kokkos::parallel_reduce( "NumSolidifiedCells", grid.domain_size_all_layers, KOKKOS_LAMBDA(const int index, int &update_meltcount, int &update_nucleatecount) { @@ -713,7 +669,7 @@ struct CellData { in_halo_region = true; if ((grain_id_all_layers_local(index) < 0) && (!in_halo_region)) update_nucleatecount++; - if ((layer_id_all_layers_local(index) != -1) && (!in_halo_region)) + if ((temperature.layer_id_all_layers(index) != -1) && (!in_halo_region)) update_meltcount++; }, melted_cells_local, nucleated_grain_cells_local); @@ -732,6 +688,32 @@ struct CellData { return vol_fraction_nucleated; } + // If cell is liquid, undercooling at solidification start is 0. If cell is active or solid, use stored value from + // temperature view + void calcUndercoolingStart(const int domain_size, Temperature &temperature) { + + auto _cell_type = cell_type; + Kokkos::parallel_for( + "calcUndercoolingStart", domain_size, KOKKOS_LAMBDA(const int &index) { + if (_cell_type(index) == Liquid) + temperature.undercooling_solidification_start(index) = 0.0; + }); + } + + // If cell is liquid or active, undercooling is 0 if superheated, otherwise calculated from the cooling + // rate/liquidus time. If cell is solid, use stored value + void calcUndercoolingCurrent(const int domain_size, Temperature &temperature, const int cycle = 0) { + + auto _cell_type = cell_type; + Kokkos::parallel_for( + "calcUndercoolingCurrent", domain_size, KOKKOS_LAMBDA(const int &index) { + if (_cell_type(index) != Solid) { + if (temperature.current_cooling_rate(index) >= 0) + temperature.undercooling_current(index) = temperature.getUndercooling(cycle, index); + } + }); + } + // If storing whether or not a cell is at a melt pool edge, update the value KOKKOS_INLINE_FUNCTION void setMeltEdge(const int index, const bool updated_state) const { @@ -743,7 +725,6 @@ struct CellData { // corresponding to the current layer of a multilayer problem auto getGrainIDSubview(const Grid &grid) const { return Kokkos::subview(grain_id_all_layers, grid.layer_range); } auto getPhaseIDSubview(const Grid &grid) const { return Kokkos::subview(phase_id_all_layers, grid.layer_range); } - auto getLayerIDSubview(const Grid &grid) const { return Kokkos::subview(layer_id_all_layers, grid.layer_range); } }; #endif diff --git a/src/CAnucleation.hpp b/src/CAnucleation.hpp index 7b727178..de39ef3a 100644 --- a/src/CAnucleation.hpp +++ b/src/CAnucleation.hpp @@ -29,6 +29,8 @@ struct Nucleation { using memory_space = MemorySpace; using view_type_int = Kokkos::View; using view_type_int_host = typename view_type_int::host_mirror_type; + using view_type_double_3d = Kokkos::View; + using view_type_double_3d_host = typename view_type_double_3d::host_mirror_type; // Using the default exec space for this memory space. using execution_space = typename memory_space::execution_space; @@ -80,21 +82,103 @@ struct Nucleation { successful_nucleation_counter = 0; } + // Get time-temperature history in the form needed for nucleation initialization from the structs stored in + // temperature + view_type_double_3d_host getTimeTempHistory(std::string simulation_type, const int layernumber, const double deltat, + const Grid &grid, const Temperature &temperature) { + view_type_double_3d_host time_temp_history_host("time_temp_history", grid.domain_size, + temperature.max_num_solidification_events, 3); + if ((simulation_type == "SingleGrain") || (simulation_type == "Directional")) { + // Simple time-temperature history, only 0-1 solidification event per cell TODO: Does not work for spot + // problem + auto time_temp_history = Kokkos::create_mirror_view(memory_space(), time_temp_history_host); + Kokkos::parallel_for( + "FillTimeTempHistory", grid.domain_size, KOKKOS_LAMBDA(const int &index) { + if (temperature.current_cooling_rate(index) > 0) { + time_temp_history(index, 0, 0) = -1.0; + time_temp_history(index, 0, 1) = temperature.last_time_below_liquidus(index); + time_temp_history(index, 0, 2) = temperature.current_cooling_rate(index); + } + }); + Kokkos::fence(); + time_temp_history_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), time_temp_history); + } + else { + // Temporary view to count events stored per cell + view_type_int_host number_of_solidification_events_host("num_solidification_events", grid.domain_size); + // First, get number of solidification events per cell and the max number of times a cell solidifies + // Data was already read into the "raw_temperature_data" data structure + // Determine which section of "raw_temperature_data" is relevant for this layer of the overall domain + int start_range = temperature.first_value[layernumber]; + int end_range = temperature.last_value[layernumber]; + for (int i = start_range; i < end_range; i++) { + // Get the integer X, Y, Z coordinates associated with this data point, along with the associated TM, + // TL, CR values coord_y is relative to this MPI rank's grid, while coord_y_global is relative to the + // overall simulation domain + int coord_x = temperature.getTempCoordX(i, grid.x_min, grid.deltax); + int coord_y = temperature.getTempCoordY(i, grid.y_min, grid.deltax, grid.y_offset); + int coord_z = + temperature.getTempCoordZ(i, grid.deltax, grid.layer_height, layernumber, grid.z_min_layer); + double t_melting = temperature.getTempCoordTM(i); + double t_liquidus = temperature.getTempCoordTL(i); + double cooling_rate_cell = temperature.getTempCoordCR(i); + + // 1D cell coordinate on this MPI rank's domain + int index = grid.get1DIndex(coord_x, coord_y, coord_z); + // Store TM, TL, CR values for this solidification event in time_temp_history_host as double values + const int n_solidification_events_cell = number_of_solidification_events_host(index); + time_temp_history_host(index, n_solidification_events_cell, 0) = + static_cast(Kokkos::round(t_melting / deltat) + 1); + time_temp_history_host(index, n_solidification_events_cell, 1) = + static_cast(Kokkos::round(t_liquidus / deltat) + 1); + // Cannot go above/below liquidus on same time step - must offset by at least 1 + if (time_temp_history_host(index, n_solidification_events_cell, 0) == + time_temp_history_host(index, n_solidification_events_cell, 1)) + time_temp_history_host(index, n_solidification_events_cell, 1)++; + time_temp_history_host(index, n_solidification_events_cell, 2) = std::abs(cooling_rate_cell) * deltat; + // Increment number of solidification events for this cell + number_of_solidification_events_host(index)++; + } + + // Reorder solidification events in time_temp_history_host(location,event number,component) so that they are + // in order based on the melting time values (component = 0 in time_temp_history_host) + for (int index = 0; index < grid.domain_size; index++) { + int n_solidification_events_cell = number_of_solidification_events_host(index); + if (n_solidification_events_cell > 0) { + for (int i = 0; i < n_solidification_events_cell - 1; i++) { + for (int j = (i + 1); j < n_solidification_events_cell; j++) { + if (time_temp_history_host(index, i, 0) > time_temp_history_host(index, j, 0)) { + // Swap these two points - melting event "j" happens before event "i" + float old_melt_val = time_temp_history_host(index, i, 0); + float old_liq_val = time_temp_history_host(index, i, 1); + float old_cr_val = time_temp_history_host(index, i, 2); + time_temp_history_host(index, i, 0) = time_temp_history_host(index, j, 0); + time_temp_history_host(index, i, 1) = time_temp_history_host(index, j, 1); + time_temp_history_host(index, i, 2) = time_temp_history_host(index, j, 2); + time_temp_history_host(index, j, 0) = old_melt_val; + time_temp_history_host(index, j, 1) = old_liq_val; + time_temp_history_host(index, j, 2) = old_cr_val; + } + } + } + } + } + } + return time_temp_history_host; + } + // Initialize nucleation site locations, GrainID values, and time at which nucleation events will potentially occur, // accounting for multiple possible nucleation events in cells that melt and solidify multiple times template - void placeNuclei(const Temperature &temperature, const InterfacialResponseFunction &irf, - const unsigned long rng_seed, const int layernumber, const Grid &grid, const int id) { + void placeNuclei(std::string simulation_type, const Temperature &temperature, const InterfacialResponseFunction &irf, const unsigned long rng_seed, const int layernumber, const Grid &grid, const int id, const double deltat) { // TODO: convert this subroutine into kokkos kernels, rather than copying data back to the host, and nucleation // data back to the device again. This is currently performed on the device due to heavy usage of standard // library algorithm functions Copy temperature data into temporary host views for this subroutine - auto max_solidification_events_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.max_solidification_events); - auto number_of_solidification_events_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.number_of_solidification_events); - auto liquidus_time_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.liquidus_time); - auto cooling_rate_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.cooling_rate); + + // Organize solidification data by event to assign nucleation to specific times that cells go below the liquidus + view_type_double_3d_host time_temp_history_host = + getTimeTempHistory(simulation_type, layernumber, deltat, grid, temperature); // Use new RNG seed for each layer std::mt19937_64 generator(rng_seed + static_cast(layernumber)); @@ -115,7 +199,7 @@ struct Nucleation { // If each cell underwent solidification 1x, the number of potential nuclei in the layer const long int nuclei_this_layer_single_long = std::lround(_inputs.n_max * domain_volume); // Multiplier for the number of nucleation events per layer, based on the max number of solidification events - long int nuclei_multiplier_long = static_cast(max_solidification_events_host(layernumber)); + long int nuclei_multiplier_long = static_cast(temperature.max_num_solidification_events); long int nuclei_this_layer_long = nuclei_this_layer_single_long * nuclei_multiplier_long; // Vector and view sizes should be type int for indexing and grain ID assignment purposes - // Nuclei_ThisLayer_long should be less than INT_MAX @@ -183,16 +267,16 @@ struct Nucleation { grid.get1DIndex(nuclei_x(n_event), nuclei_y(n_event) - grid.y_offset, nuclei_z(n_event)); // Criteria for placing a nucleus - whether or not this nuclei is associated with a solidification // event - if (meltevent < number_of_solidification_events_host(nuclei_location_this_layer)) { + if (meltevent < temperature.number_of_solidification_events(nuclei_location_this_layer)) { // Nucleation event is possible - cell undergoes solidification at least once, this nucleation // event is associated with one of the time periods during which the associated cell undergoes // solidification nuclei_location_myrank_v[possible_nuclei] = nuclei_location_this_layer; double liq_time_this_event = - static_cast(liquidus_time_host(nuclei_location_this_layer, meltevent, 1)); + static_cast(time_temp_history_host(nuclei_location_this_layer, meltevent, 1)); double cooling_rate_this_event = - static_cast(cooling_rate_host(nuclei_location_this_layer, meltevent)); + static_cast(time_temp_history_host(nuclei_location_this_layer, meltevent, 2)); double time_to_nuc_und = liq_time_this_event + nuclei_undercooling_whole_domain_v[n_event] / cooling_rate_this_event; if (liq_time_this_event > time_to_nuc_und) diff --git a/src/CAprint.hpp b/src/CAprint.hpp index ab6eeec6..b200f099 100644 --- a/src/CAprint.hpp +++ b/src/CAprint.hpp @@ -206,8 +206,8 @@ struct Print { const Grid &grid, CellData &celldata, Temperature &temperature, Interface &interface, Orientation &orientation) { - using view_type_float = Kokkos::View; - using view_type_int = Kokkos::View; + using view_type_float_host = Kokkos::View; + using view_type_int_host = Kokkos::View; if ((_inputs.intralayer) && (cycle % _inputs.intralayer_increment == 0)) { // Current time in microseconds intralayer_times.push_back(getCurrentTime(cycle, deltat)); @@ -226,8 +226,7 @@ struct Print { printViewData(id, intralayer_ofstream, grid, true, "int", "GrainID", grain_id_whole_domain); } if (_inputs.intralayer_layer_id) { - auto layer_id_current_layer = celldata.getLayerIDSubview(grid); - auto layer_id_whole_domain = collectViewData(id, np, grid, true, MPI_SHORT, layer_id_current_layer); + auto layer_id_whole_domain = collectViewData(id, np, grid, true, MPI_SHORT, temperature.layer_id); printViewData(id, intralayer_ofstream, grid, true, "short", "LayerID", layer_id_whole_domain); } if (_inputs.intralayer_phase_id) { @@ -236,32 +235,33 @@ struct Print { printViewData(id, intralayer_ofstream, grid, true, "short", "PhaseID", phase_id_whole_domain); } if (_inputs.intralayer_undercooling_current) { - // TODO: remove this, the subview undercooling_current_layer is already stored in the temperature struct - auto undercooling_current_layer = - Kokkos::subview(temperature.undercooling_current_all_layers, grid.layer_range); + // Get undercooling data, accounting for cell types + celldata.calcUndercoolingCurrent(grid.domain_size, temperature, cycle); auto undercooling_whole_domain = - collectViewData(id, np, grid, true, MPI_FLOAT, undercooling_current_layer); + collectViewData(id, np, grid, true, MPI_FLOAT, temperature.undercooling_current); printViewData(id, intralayer_ofstream, grid, true, "float", "UndercoolingCurrent", undercooling_whole_domain); } if (_inputs.intralayer_undercooling_solidification_start) { + // Get undercooling data, accounting for cell types + celldata.calcUndercoolingStart(grid.domain_size, temperature); auto undercooling_start_whole_domain = collectViewData(id, np, grid, true, MPI_FLOAT, temperature.undercooling_solidification_start); printViewData(id, intralayer_ofstream, grid, true, "float", "UndercoolingStart", undercooling_start_whole_domain); } if (_inputs.intralayer_melt_time_step) { - auto melt_time_step = temperature.template extractTmTlData(0, grid.domain_size); + auto melt_time_step = temperature.template extractTmData(grid.domain_size); auto melt_time_step_whole_domain = collectViewData(id, np, grid, true, MPI_INT, melt_time_step); printViewData(id, intralayer_ofstream, grid, true, "int", "MeltTimeStep", melt_time_step_whole_domain); } if (_inputs.intralayer_crit_time_step) { - auto crit_time_step = temperature.template extractTmTlData(1, grid.domain_size); + auto crit_time_step = temperature.template extractTlData(grid.domain_size); auto crit_time_step_whole_domain = collectViewData(id, np, grid, true, MPI_INT, crit_time_step); printViewData(id, intralayer_ofstream, grid, true, "int", "CritTimeStep", crit_time_step_whole_domain); } if (_inputs.intralayer_undercooling_change) { - auto cooling_rate = temperature.template extractCrData(grid.domain_size); + auto cooling_rate = temperature.template extractCrData(grid.domain_size); auto undercooling_change_whole_domain = collectViewData(id, np, grid, true, MPI_FLOAT, cooling_rate); printViewData(id, intralayer_ofstream, grid, true, "float", "UndercoolingChange", undercooling_change_whole_domain); @@ -334,8 +334,8 @@ struct Print { if (id == 0) std::cout << "Layer " << layernumber << " finished solidification" << std::endl; - using view_type_float = Kokkos::View; - using view_type_int = Kokkos::View; + using view_type_float_host = Kokkos::View; + using view_type_int_host = Kokkos::View; if ((!_inputs.skip_all_printing) && (layernumber == _inputs.print_layer_number[interlayer_file_count])) { // If printing a time series, include this file and time as the final one if (_inputs.intralayer) { @@ -365,7 +365,7 @@ struct Print { grain_id_all_layers_whole_domain); if (_inputs.interlayer_layer_id) { auto layer_id_all_layers_whole_domain = - collectViewData(id, np, grid, false, MPI_SHORT, celldata.layer_id_all_layers); + collectViewData(id, np, grid, false, MPI_SHORT, temperature.layer_id_all_layers); printViewData(id, interlayer_all_layers_ofstream, grid, false, "short", "LayerID", layer_id_all_layers_whole_domain); } @@ -376,12 +376,16 @@ struct Print { phase_id_all_layers_whole_domain); } if (_inputs.interlayer_undercooling_current) { + // Final state of undercooling field in all cells + celldata.calcUndercoolingCurrent(grid.domain_size, temperature); auto undercooling_all_layers_whole_domain = collectViewData(id, np, grid, false, MPI_FLOAT, temperature.undercooling_current_all_layers); printViewData(id, interlayer_all_layers_ofstream, grid, false, "float", "UndercoolingCurrent", undercooling_all_layers_whole_domain); } if (_inputs.interlayer_undercooling_solidification_start) { + // Undercooling for the last time each cell underwent solidification + celldata.calcUndercoolingStart(grid.domain_size, temperature); auto undercooling_start_all_layers_whole_domain = collectViewData( id, np, grid, false, MPI_FLOAT, temperature.undercooling_solidification_start_all_layers); printViewData(id, interlayer_all_layers_ofstream, grid, false, "float", "UndercoolingStart", @@ -414,19 +418,19 @@ struct Print { writeHeader(currentlayer_ofstream, vtk_filename_current_layer, grid, true); } if (_inputs.interlayer_melt_time_step) { - auto melt_time_step = temperature.template extractTmTlData(0, grid.domain_size); + auto melt_time_step = temperature.template extractTmData(grid.domain_size); auto melt_time_step_whole_domain = collectViewData(id, np, grid, true, MPI_INT, melt_time_step); printViewData(id, currentlayer_ofstream, grid, true, "int", "MeltTimeStep", melt_time_step_whole_domain); } if (_inputs.interlayer_crit_time_step) { - auto crit_time_step = temperature.template extractTmTlData(1, grid.domain_size); + auto crit_time_step = temperature.template extractTlData(grid.domain_size); auto crit_time_step_whole_domain = collectViewData(id, np, grid, true, MPI_INT, crit_time_step); printViewData(id, currentlayer_ofstream, grid, true, "int", "CritTimeStep", crit_time_step_whole_domain); } if (_inputs.interlayer_undercooling_change) { - auto cooling_rate = temperature.template extractCrData(grid.domain_size); + auto cooling_rate = temperature.template extractCrData(grid.domain_size); auto undercooling_change_whole_domain = collectViewData(id, np, grid, true, MPI_FLOAT, cooling_rate); printViewData(id, currentlayer_ofstream, grid, true, "float", "UndercoolingChange", @@ -475,7 +479,7 @@ struct Print { // If necessary, print RVE data for ExaConstit input if (_inputs.print_default_rve) { auto layer_id_all_layers_whole_domain = - collectViewData(id, np, grid, false, MPI_SHORT, celldata.layer_id_all_layers); + collectViewData(id, np, grid, false, MPI_SHORT, temperature.layer_id_all_layers); if (id == 0) printExaConstitDefaultRVE(grid, layer_id_all_layers_whole_domain, grain_id_all_layers_whole_domain); } diff --git a/src/CAtemperature.hpp b/src/CAtemperature.hpp index bcaa6a74..d3e6bcb3 100644 --- a/src/CAtemperature.hpp +++ b/src/CAtemperature.hpp @@ -24,13 +24,17 @@ template struct Temperature { using memory_space = MemorySpace; + using view_type_short = Kokkos::View; using view_type_int = Kokkos::View; + using view_type_int_2d = Kokkos::View; using view_type_int_3d = Kokkos::View; using view_type_double = Kokkos::View; using view_type_double_2d = Kokkos::View; using view_type_float = Kokkos::View; using view_type_float_2d = Kokkos::View; + using view_type_short_host = typename view_type_short::host_mirror_type; using view_type_int_host = typename view_type_int::host_mirror_type; + using view_type_int_2d_host = typename view_type_int_2d::host_mirror_type; using view_type_int_3d_host = typename view_type_int_3d::host_mirror_type; using view_type_double_host = typename view_type_double::host_mirror_type; using view_type_double_2d_host = typename view_type_double_2d::host_mirror_type; @@ -41,20 +45,24 @@ struct Temperature { // Using the default exec space for this memory space. using execution_space = typename memory_space::execution_space; - // Maximum number of times each CA cell in a given layer undergoes solidification - view_type_int max_solidification_events; - // Each time that a cell (index 0) goes above and below the liquidus time (index 2) for each melt-solidification - // event (index 1) - view_type_int_3d liquidus_time; - // Cooling rate for each cell (index 1) for each solidification event (index 2) - view_type_float_2d cooling_rate; + int liquidus_time_counter = 0; + int num_liquidus_times_this_layer = 0; + // Only nonzero in the special case of 0 cooling rate/uniform undercooling + float init_undercooling = 0.0; + // Times at which each cell crosses the liquidus temperature + view_type_int_host liquidus_time_list_host; + // Cells associated with each value in liquidus_time_list_host + view_type_int liquidus_cell_list; + // Cooling rate associated with each value in liquidus_time_list_host (INT_MAX placeholder for cells that are + // heating) + view_type_float cooling_rate_list; + // The last time a given cell cooled below the liquidus (0.0 placeholder if it hasn't cooled below the liquidus) + view_type_int last_time_below_liquidus; + // Cooling rate currently in use by each cell + view_type_float current_cooling_rate; // The number of times that each CA cell will traverse the liquidus during this layer - view_type_int number_of_solidification_events; - // A counter for the number of times each CA cell has undergone solidification so far this layer - view_type_int solidification_event_counter; - // The current undercooling of a CA cell (if superheated liquid or hasn't undergone solidification yet, equals 0) - // Also maintained for the full multilayer domain - view_type_float undercooling_current_all_layers, undercooling_current; + view_type_int_host number_of_solidification_events; + int max_num_solidification_events = 1; // Data structure for storing raw temperature data from file(s) // Store data as double - needed for small time steps to resolve local differences in solidification conditions // Each data point has 6 values (X, Y, Z coordinates, melting time, liquidus time, and cooling rate) @@ -63,70 +71,116 @@ struct Temperature { // These contain "number_of_layers" values corresponding to the location within "raw_temperature_data" of the first // data element in each temperature file, if used view_type_int_host first_value, last_value; + // Layer associated with last time a cell undergoes melting/soldification + view_type_short layer_id_all_layers, layer_id; + // Number of times in the layer a cell has undergone solidification so far (optional to store based on selected + // inputs) + bool store_solidification_event_counter = false; + view_type_int solidification_event_counter; // Undercooling when a solidification first began in a cell (optional to store based on selected inputs) - bool _store_solidification_start; + bool store_undercooling_current = false; + view_type_float undercooling_current_all_layers, undercooling_current; + // Undercooling when a solidification ended in a cell (optional to store based on selected inputs) + bool store_undercooling_solidification_start = false; view_type_float undercooling_solidification_start_all_layers, undercooling_solidification_start; // Temperature field inputs from file TemperatureInputs _inputs; - // Constructor creates views with size based on the grid inputs - liquidus_time and cooling_rate are later modified - // to account for multiple events if needed, undercooling_current and solidification_event_counter are default - // initialized to zeros - Temperature(const Grid &grid, TemperatureInputs inputs, const bool store_solidification_start = false, + // Constructor creates views with size based on the grid inputs. If any of the output options + // solidification_event_counter, undercooling_solidification_start, undercooling_current, and layer_id are toggled, + // the associated views are allocated and filled during the simulation + Temperature(const Grid &grid, TemperatureInputs inputs, PrintInputs print_inputs, const int est_num_temperature_data_points = 100000) - : max_solidification_events( - view_type_int(Kokkos::ViewAllocateWithoutInitializing("number_of_layers"), grid.number_of_layers)) - , liquidus_time( - view_type_int_3d(Kokkos::ViewAllocateWithoutInitializing("liquidus_time"), grid.domain_size, 1, 2)) - , cooling_rate(view_type_float_2d(Kokkos::ViewAllocateWithoutInitializing("cooling_rate"), grid.domain_size, 1)) - , number_of_solidification_events(view_type_int( + : liquidus_time_list_host( + view_type_int_host(Kokkos::ViewAllocateWithoutInitializing("liquidus_time_list"), grid.domain_size)) + , liquidus_cell_list( + view_type_int(Kokkos::ViewAllocateWithoutInitializing("liquidus_cell_list"), grid.domain_size)) + , cooling_rate_list( + view_type_float(Kokkos::ViewAllocateWithoutInitializing("cooling_rate_list"), grid.domain_size)) + , last_time_below_liquidus(view_type_int("last_time_below_liquidus", grid.domain_size)) + , current_cooling_rate(view_type_float("current_cooling_rate", grid.domain_size)) + , number_of_solidification_events(view_type_int_host( Kokkos::ViewAllocateWithoutInitializing("number_of_solidification_events"), grid.domain_size)) - , solidification_event_counter(view_type_int("solidification_event_counter", grid.domain_size)) - , undercooling_current_all_layers(view_type_float("undercooling_current", grid.domain_size_all_layers)) , raw_temperature_data(view_type_double_2d_host(Kokkos::ViewAllocateWithoutInitializing("raw_temperature_data"), est_num_temperature_data_points, num_temperature_components)) , first_value(view_type_int_host(Kokkos::ViewAllocateWithoutInitializing("first_value"), grid.number_of_layers)) , last_value(view_type_int_host(Kokkos::ViewAllocateWithoutInitializing("last_value"), grid.number_of_layers)) - , _store_solidification_start(store_solidification_start) + , layer_id_all_layers(view_type_short(Kokkos::ViewAllocateWithoutInitializing("layer_id_all_layers"), + grid.domain_size_all_layers)) + , solidification_event_counter( + view_type_int(Kokkos::ViewAllocateWithoutInitializing("solidification_event_counter"), 1)) + , undercooling_current_all_layers(view_type_float( + Kokkos::ViewAllocateWithoutInitializing("undercooling_current_all_layers"), grid.domain_size_all_layers)) + , undercooling_solidification_start_all_layers( + view_type_float(Kokkos::ViewAllocateWithoutInitializing("undercooling_current_all_layers"), 1)) + , _inputs(inputs) { - if (_store_solidification_start) { - // Default init starting undercooling in cells to zero - undercooling_solidification_start_all_layers = - view_type_float("undercooling_solidification_start", grid.domain_size_all_layers); - getCurrentLayerStartingUndercooling(grid.layer_range); - } - getCurrentLayerUndercooling(grid.layer_range); + Kokkos::deep_copy(layer_id_all_layers, -1); + layer_id = getLayerSubview(layer_id_all_layers, grid.layer_range); + + // Allocate views for optional stored data based on print inputs + initOptionalViews(grid, print_inputs); } // Constructor using in-memory temperature data from external source (input_temperature_data) - Temperature(const int id, const int np, const Grid &grid, TemperatureInputs inputs, - view_type_coupled input_temperature_data, const bool store_solidification_start = false) - : max_solidification_events( - view_type_int(Kokkos::ViewAllocateWithoutInitializing("number_of_layers"), grid.number_of_layers)) - , liquidus_time( - view_type_int_3d(Kokkos::ViewAllocateWithoutInitializing("liquidus_time"), grid.domain_size, 1, 2)) - , cooling_rate(view_type_float_2d(Kokkos::ViewAllocateWithoutInitializing("cooling_rate"), grid.domain_size, 1)) - , number_of_solidification_events(view_type_int( + Temperature(const int id, const int np, const Grid &grid, TemperatureInputs inputs, PrintInputs print_inputs, + view_type_coupled input_temperature_data, const int est_num_temperature_data_points = 100000) + : liquidus_time_list_host( + view_type_int_host(Kokkos::ViewAllocateWithoutInitializing("liquidus_time_list"), 2 * grid.domain_size)) + , liquidus_cell_list( + view_type_int(Kokkos::ViewAllocateWithoutInitializing("liquidus_cell_list"), 2 * grid.domain_size)) + , cooling_rate_list( + view_type_float(Kokkos::ViewAllocateWithoutInitializing("cooling_rate_list"), 2 * grid.domain_size)) + , last_time_below_liquidus(view_type_int("last_time_below_liquidus", grid.domain_size)) + , current_cooling_rate(view_type_float("current_cooling_rate", grid.domain_size)) + , number_of_solidification_events(view_type_int_host( Kokkos::ViewAllocateWithoutInitializing("number_of_solidification_events"), grid.domain_size)) - , solidification_event_counter(view_type_int("solidification_event_counter", grid.domain_size)) - , undercooling_current_all_layers(view_type_float("undercooling_current", grid.domain_size_all_layers)) , raw_temperature_data(view_type_double_2d_host(Kokkos::ViewAllocateWithoutInitializing("raw_temperature_data"), - input_temperature_data.extent(0) * (inputs.number_of_copies), - num_temperature_components)) + est_num_temperature_data_points, num_temperature_components)) , first_value(view_type_int_host(Kokkos::ViewAllocateWithoutInitializing("first_value"), grid.number_of_layers)) , last_value(view_type_int_host(Kokkos::ViewAllocateWithoutInitializing("last_value"), grid.number_of_layers)) - , _store_solidification_start(store_solidification_start) + , layer_id_all_layers(view_type_short(Kokkos::ViewAllocateWithoutInitializing("layer_id_all_layers"), + grid.domain_size_all_layers)) + , solidification_event_counter( + view_type_int(Kokkos::ViewAllocateWithoutInitializing("solidification_event_counter"), 1)) + , undercooling_current_all_layers(view_type_float( + Kokkos::ViewAllocateWithoutInitializing("undercooling_current_all_layers"), grid.domain_size_all_layers)) + , undercooling_solidification_start_all_layers( + view_type_float(Kokkos::ViewAllocateWithoutInitializing("undercooling_current_all_layers"), 1)) , _inputs(inputs) { copyTemperatureData(id, np, grid, input_temperature_data); - if (_store_solidification_start) { - // Default init starting undercooling in cells to zero - undercooling_solidification_start_all_layers = - view_type_float("undercooling_solidification_start", grid.domain_size_all_layers); - getCurrentLayerStartingUndercooling(grid.layer_range); + initOptionalViews(grid, print_inputs); + Kokkos::deep_copy(layer_id_all_layers, -1); + layer_id = getLayerSubview(layer_id_all_layers, grid.layer_range); + } + + void initOptionalViews(const Grid &grid, PrintInputs print_inputs) { + if ((print_inputs.intralayer_solidification_event_counter) || + (print_inputs.interlayer_solidification_event_counter)) { + store_solidification_event_counter = true; + // Solidification event counter is only stored for the current layer, init to zeros + Kokkos::realloc(solidification_event_counter, grid.domain_size); + Kokkos::deep_copy(solidification_event_counter, 0); + } + if ((print_inputs.intralayer_undercooling_current) || (print_inputs.interlayer_undercooling_current)) { + store_undercooling_current = true; + // Undercooling is stored for all layers, as is the subview for the current layer's undercooling, init both + // to zeros + Kokkos::realloc(undercooling_current_all_layers, grid.domain_size_all_layers); + Kokkos::deep_copy(undercooling_current_all_layers, 0.0); + // Subview for current layer + undercooling_current = getLayerSubview(undercooling_current_all_layers, grid.layer_range); + } + if ((print_inputs.intralayer_undercooling_solidification_start) || + (print_inputs.interlayer_undercooling_solidification_start)) { + store_undercooling_solidification_start = true; + Kokkos::realloc(undercooling_solidification_start_all_layers, grid.domain_size_all_layers); + Kokkos::deep_copy(undercooling_solidification_start_all_layers, 0.0); + undercooling_solidification_start = + getLayerSubview(undercooling_solidification_start_all_layers, grid.layer_range); } - getCurrentLayerUndercooling(grid.layer_range); } // If using a problem type that involves translating or mirroring temperature data, output a new X or Y value based @@ -443,60 +497,105 @@ struct Temperature { else location_liquidus_isotherm = location_init_undercooling + Kokkos::round(_inputs.init_undercooling / (_inputs.G * grid.deltax)); - - // Local copies for lambda capture. - auto liquidus_time_local = liquidus_time; - auto cooling_rate_local = cooling_rate; - auto max_solidification_events_local = max_solidification_events; - auto number_solidification_events_local = number_of_solidification_events; - auto undercooling_current_local = undercooling_current; - double init_undercooling_local = _inputs.init_undercooling; - double G_local = _inputs.G; - double R_local = _inputs.R; - auto policy = Kokkos::RangePolicy(0, grid.domain_size); - Kokkos::parallel_for( - "TempInitG", policy, KOKKOS_LAMBDA(const int &index) { - // All cells past melting time step - liquidus_time_local(index, 0, 0) = -1; - // Negative dist_from_liquidus and dist_from_init_undercooling values for cells below the liquidus - // isotherm - int coord_z = grid.getCoordZ(index); - int dist_from_liquidus = coord_z - location_liquidus_isotherm; - // Cells reach liquidus at a time dependent on their Z coordinate - // Cells with negative liquidus time values are already undercooled, should have positive undercooling - // and negative liquidus time step - if (dist_from_liquidus < 0) { - liquidus_time_local(index, 0, 1) = -1; - int dist_from_init_undercooling = coord_z - location_init_undercooling; - undercooling_current_local(index) = - init_undercooling_local - dist_from_init_undercooling * G_local * grid.deltax; - } - else { - // Cells with positive liquidus time values are not yet tracked - leave current undercooling as - // default zeros and set liquidus time step. R_local will never be zero here, as all cells at a - // fixed undercooling must be below the liquidus (distFromLiquidus < 0) - liquidus_time_local(index, 0, 1) = dist_from_liquidus * G_local * grid.deltax / (R_local * deltat); + // If cooling rate is 0, need to set an initial undercooling + if (_inputs.R == 0) + init_undercooling = _inputs.init_undercooling; + + // List of liquidus time events - Since iteration over liqudius time events does not occur for the SingleGrain + // and Directional problem types (as they have a single last_time_below_liquidus and current_cooling_rate that + // apply to all cells, and the Active cells present are undercooled at the start of the simulation), the + // liquidus_time_list_host, liquidus_cell_list, and cooling_rate_list views are filled with data only for + // reference in printing temperature field info + std::vector liquidus_time_step_read(2 * grid.domain_size), cell_read(2 * grid.domain_size); + std::vector cooling_rate_read(2 * grid.domain_size); + auto current_cooling_rate_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), current_cooling_rate); + auto last_time_below_liquidus_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), last_time_below_liquidus); + for (int coord_z = 0; coord_z < grid.nz; coord_z++) { + const int dist_from_liquidus = coord_z - location_liquidus_isotherm; + for (int coord_x = 0; coord_x < grid.nx; coord_x++) { + for (int coord_y = 0; coord_y < grid.ny_local; coord_y++) { + // 1D cell index + const int index = grid.get1DIndex(coord_x, coord_y, coord_z); + // Initially melted (time step 0) + liquidus_time_step_read[num_liquidus_times_this_layer] = 0; + cell_read[num_liquidus_times_this_layer] = index; + cooling_rate_read[num_liquidus_times_this_layer] = -1.0; + num_liquidus_times_this_layer++; + // Cells reach liquidus at a time dependent on their Z coordinate + // Cells with negative last_time_below_liquidus are already undercooled, times may be negative + liquidus_time_step_read[num_liquidus_times_this_layer] = + Kokkos::round(dist_from_liquidus * _inputs.G * grid.deltax / (_inputs.R * deltat)); + cell_read[num_liquidus_times_this_layer] = index; + // Cells cool at a constant rate + cooling_rate_read[num_liquidus_times_this_layer] = _inputs.R * deltat; + // Cell will only have one cooling rate/liquidus time + last_time_below_liquidus_host(index) = liquidus_time_step_read[num_liquidus_times_this_layer]; + current_cooling_rate_host(index) = cooling_rate_read[num_liquidus_times_this_layer]; + num_liquidus_times_this_layer++; } - // Cells cool at a constant rate - cooling_rate_local(index, 0) = R_local * deltat; - // All cells solidify once - max_solidification_events_local(0) = 1; - number_solidification_events_local(index) = 1; - }); + } + } + + // From the 3 vectors liquidus_time_step_read, cell_read, cooling_rate_read initialize the views + // liquidus_time_list_host, liquidus_cell_list (stored on device), and cooling_rate_list (stored on device) + initOrderedTimeTempHistory(liquidus_time_step_read, cell_read, cooling_rate_read, + num_liquidus_times_this_layer); + + // Counter goes unused since all melting events have occurred and the time at which the cells go below the + // liquidus/their cooling rates are also known. Set counter to max value + liquidus_time_counter = num_liquidus_times_this_layer; + + // Copy host view data back to device + current_cooling_rate = Kokkos::create_mirror_view_and_copy(memory_space(), current_cooling_rate_host); + last_time_below_liquidus = Kokkos::create_mirror_view_and_copy(memory_space(), last_time_below_liquidus_host); + + // Each cell solidifies once + Kokkos::deep_copy(number_of_solidification_events, 1); + Kokkos::deep_copy(layer_id_all_layers, 0); MPI_Barrier(MPI_COMM_WORLD); if (id == 0) { - std::cout << "Temperature field initialized for unidirectional solidification with G = " << G_local + std::cout << "Temperature field initialized for unidirectional solidification with G = " << _inputs.G << " K/m, initial undercooling at Z = " << location_init_undercooling << " of " << _inputs.init_undercooling << " K below the liquidus" << std::endl; std::cout << "Done with temperature field initialization" << std::endl; } } - // Initialize temperature data for a hemispherical spot melt (single layer simulation) - void initialize(const int id, const Grid &grid, const double freezing_range, double deltat, double spot_radius) { + // From the 3 vectors liquidus_time_step_vec, cell_vec, cooling_rate_vec, initialize the views + // liquidus_time_list_host, liquidus_cell_list (stored on device), and cooling_rate_list (stored on device) + void initOrderedTimeTempHistory(std::vector liquidus_time_step_vec, std::vector cell_vec, + std::vector cooling_rate_vec, const int number_vals_stored) { + + // Reorder liquidus time events based on the times at which they occur + std::vector> time_temp_history; + time_temp_history.reserve(number_vals_stored); + for (int n = 0; n < number_vals_stored; n++) { + time_temp_history.push_back(std::make_tuple(liquidus_time_step_vec[n], cell_vec[n], cooling_rate_vec[n])); + } + std::sort(time_temp_history.begin(), time_temp_history.end()); + + // Resize views for number of values in the vectors + Kokkos::realloc(liquidus_time_list_host, number_vals_stored); + Kokkos::realloc(liquidus_cell_list, number_vals_stored); + Kokkos::realloc(cooling_rate_list, number_vals_stored); + + // Copy to host view for liquidus_time_list, temporary host views for list of cells and cooling rates + auto liquidus_cell_list_host = Kokkos::create_mirror_view(Kokkos::HostSpace(), liquidus_cell_list); + auto cooling_rate_list_host = Kokkos::create_mirror_view(Kokkos::HostSpace(), cooling_rate_list); + for (int n = 0; n < num_liquidus_times_this_layer; n++) { + liquidus_time_list_host(n) = std::get<0>(time_temp_history[n]); + liquidus_cell_list_host(n) = std::get<1>(time_temp_history[n]); + cooling_rate_list_host(n) = std::get<2>(time_temp_history[n]); + } + + // Copy host view data back to device + liquidus_cell_list = Kokkos::create_mirror_view_and_copy(memory_space(), liquidus_cell_list_host); + cooling_rate_list = Kokkos::create_mirror_view_and_copy(memory_space(), cooling_rate_list_host); + } - // Each cell solidifies at most one time - Kokkos::deep_copy(max_solidification_events, 1); + // Store temperature data for a hemispherical spot melt (single layer simulation) + void storeSpotData(const int id, const Grid &grid, const double freezing_range, double deltat, double spot_radius) { // Outer edges of spots are initialized at the liquidus temperature // Spots cool at constant rate R, spot thermal gradient = G @@ -511,171 +610,145 @@ struct Temperature { std::cout << "Initializing temperature field for the hemispherical spot, which will take approximately " << spot_time_est << " time steps to fully solidify" << std::endl; - // Initialize liquidus_time and cooling_rate values for this spot - auto liquidus_time_local = liquidus_time; - auto cooling_rate_local = cooling_rate; - auto number_of_solidification_events_local = number_of_solidification_events; - double R_local = _inputs.R; - auto md_policy = - Kokkos::MDRangePolicy>( - {0, 0, 0}, {grid.nz, grid.nx, grid.ny_local}); - Kokkos::parallel_for( - "SpotTemperatureInit", md_policy, KOKKOS_LAMBDA(const int coord_z, const int coord_x, const int coord_y) { - // 1D cell index - int index = grid.get1DIndex(coord_x, coord_y, coord_z); - // Distance of this cell from the spot center - float dist_z = spot_center_z - coord_z; - float dist_x = spot_center_x - coord_x; - int coord_y_global = coord_y + grid.y_offset; - float dist_y = spot_center_y - coord_y_global; - float tot_dist = Kokkos::hypot(dist_x, dist_y, dist_z); - if (tot_dist <= spot_radius) { - // Melt time step (first time step) - liquidus_time_local(index, 0, 0) = 1; - // Liquidus time step (related to distance to spot edge). Must be at least 1 time step after melt - // time step - liquidus_time_local(index, 0, 1) = Kokkos::round((spot_radius - tot_dist) / isotherm_velocity) + 2; - // Cooling rate per time step - cooling_rate_local(index, 0) = R_local * deltat; - number_of_solidification_events_local(index) = 1; - } - else { - // Dummy layer_time_temp_history data - liquidus_time_local(index, 0, 0) = -1; - liquidus_time_local(index, 0, 1) = -1; - cooling_rate_local(index, 0) = -1.0; - number_of_solidification_events_local(index) = 0; + // Add spot data to raw_temperature_data as if read from a file + // Ensure raw_temperature_data can store all data + Kokkos::realloc(raw_temperature_data, grid.domain_size, 6); + // Iterate over cells on this rank, determine if inside of spot + for (int coord_z = 0; coord_z < grid.nz; coord_z++) { + for (int coord_x = 0; coord_x < grid.nx; coord_x++) { + for (int coord_y = 0; coord_y < grid.ny_local; coord_y++) { + // Distance of this cell from the spot center + float dist_z = spot_center_z - coord_z; + float dist_x = spot_center_x - coord_x; + int coord_y_global = coord_y + grid.y_offset; + float dist_y = spot_center_y - coord_y_global; + float tot_dist = Kokkos::hypot(dist_x, dist_y, dist_z); + if (tot_dist <= spot_radius) { + // x, y, z of location in global domain + raw_temperature_data(num_liquidus_times_this_layer, 0) = + grid.x_min + static_cast(coord_x) * grid.deltax; + raw_temperature_data(num_liquidus_times_this_layer, 1) = + grid.y_min + static_cast(coord_y_global) * grid.deltax; + raw_temperature_data(num_liquidus_times_this_layer, 2) = + grid.z_min + static_cast(coord_z) * grid.deltax; + // Melting time in seconds - all cells melt at the start of the simulation + raw_temperature_data(num_liquidus_times_this_layer, 3) = 0.0; + // Liquidus time in seconds (related to distance to spot edge), must be >= deltat + raw_temperature_data(num_liquidus_times_this_layer, 4) = + deltat * (((spot_radius - tot_dist) / isotherm_velocity) + 1.0); + // Cooling rate per time step (K/s) + raw_temperature_data(num_liquidus_times_this_layer, 5) = _inputs.R; + // Increment counter for points on this rank that melt/resolidify + num_liquidus_times_this_layer++; + } } - }); + } + } + // Remove empty space in raw_temperature_data + Kokkos::resize(raw_temperature_data, num_liquidus_times_this_layer, 6); + // Start/end indices of temperature data in raw_temperature_data + first_value(0) = 0; + last_value(0) = num_liquidus_times_this_layer; MPI_Barrier(MPI_COMM_WORLD); if (id == 0) - std::cout << "Spot melt temperature field initialized" << std::endl; + std::cout + << "Spot melt temperature field initialized, number of cells to undergo melting/solidification is " + << num_liquidus_times_this_layer << std::endl; } // Calculate the number of times that a cell in layer "layernumber" undergoes melting/solidification, and store in // max_solidification_events_host - void calcMaxSolidificationEvents(const int id, const int layernumber, - const view_type_int_host max_solidification_events_host, const int start_range, - const int end_range, const Grid &grid, const std::string simulation_type) { - - bool calc_remelting_events; - if (simulation_type == "FromFile") { - if (layernumber > _inputs.temp_files_in_series) - calc_remelting_events = false; - else - calc_remelting_events = true; - } - else { - if (layernumber > 0) - calc_remelting_events = false; - else - calc_remelting_events = true; - } - if (!calc_remelting_events) { - // Use the value from a previously checked layer, since the time-temperature history is reused - if ((simulation_type == "FromFinch") || (_inputs.temp_files_in_series == 1)) { - // All layers have the same temperature data, max_solidification_events for this layer is the same as - // the last - max_solidification_events_host(layernumber) = max_solidification_events_host(layernumber - 1); - } - else { - // All layers have different temperature data but in a repeating pattern - int repeated_file = layernumber % _inputs.temp_files_in_series; - max_solidification_events_host(layernumber) = max_solidification_events_host(repeated_file); - } + view_type_int_host calcNumberOfSolidificationEvents(const int id, const int layernumber, const int start_range, + const int end_range, const Grid &grid) { + + view_type_int_host _number_of_solidification_events("number_of_solidification_events_temp", grid.domain_size); + for (int i = start_range; i < end_range; i++) { + // Get the integer X, Y, Z coordinates associated with this data point, with the Y coordinate based on + // local MPI rank's grid + int coord_x = getTempCoordX(i, grid.x_min, grid.deltax); + int coord_y = getTempCoordY(i, grid.y_min, grid.deltax, grid.y_offset); + int coord_z = getTempCoordZ(i, grid.deltax, grid.layer_height, layernumber, grid.z_min_layer); + // Convert to 1D coordinate in the current layer's domain + int index = grid.get1DIndex(coord_x, coord_y, coord_z); + _number_of_solidification_events(index)++; } - else { - // Need to calculate max_solidification_events(layernumber) - // Init to 0 - view_type_int_host temp_melt_count("temp_melt_count", grid.domain_size); - - for (int i = start_range; i < end_range; i++) { - - // Get the integer X, Y, Z coordinates associated with this data point, with the Y coordinate based on - // local MPI rank's grid - int coord_x = getTempCoordX(i, grid.x_min, grid.deltax); - int coord_y = getTempCoordY(i, grid.y_min, grid.deltax, grid.y_offset); - int coord_z = getTempCoordZ(i, grid.deltax, grid.layer_height, layernumber, grid.z_min_layer); - // Convert to 1D coordinate in the current layer's domain - int index = grid.get1DIndex(coord_x, coord_y, coord_z); - temp_melt_count(index)++; - } - int max_count = 0; - for (int i = 0; i < grid.domain_size; i++) { - if (temp_melt_count(i) > max_count) - max_count = temp_melt_count(i); - } - int max_count_global; - MPI_Allreduce(&max_count, &max_count_global, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); - max_solidification_events_host(layernumber) = max_count_global; + // Get the max of _number_of_solidification_events across all ranks + int max_count = 0; + for (int i = 0; i < grid.domain_size; i++) { + if (_number_of_solidification_events(i) > max_count) + max_count = _number_of_solidification_events(i); } - MPI_Barrier(MPI_COMM_WORLD); + MPI_Allreduce(&max_count, &max_num_solidification_events, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD); if (id == 0) std::cout << "The maximum number of melting/solidification events during layer " << layernumber << " is " - << max_solidification_events_host(layernumber) << std::endl; + << max_num_solidification_events << std::endl; + return _number_of_solidification_events; } // Read data from storage, and calculate the normalized x value of the data point - int getTempCoordX(const int i, const double x_min, const double deltax) { + int getTempCoordX(const int i, const double x_min, const double deltax) const { int x_coord = Kokkos::round((raw_temperature_data(i, 0) - x_min) / deltax); return x_coord; } // Read data from storage, and calculate the normalized y value of the data point. If the optional offset argument // is given, the return value is calculated relative to the edge of the MPI rank's local simulation domain (which is // offset by y_offset cells from the global domain edge) - int getTempCoordY(const int i, const double y_min, const double deltax, const int y_offset = 0) { + int getTempCoordY(const int i, const double y_min, const double deltax, const int y_offset = 0) const { int y_coord = Kokkos::round((raw_temperature_data(i, 1) - y_min) / deltax) - y_offset; return y_coord; } // Read data from storage, and calculate the normalized z value of the data point int getTempCoordZ(const int i, const double deltax, const int layer_height, const int layer_counter, - const view_type_double_host z_min_layer) { + const view_type_double_host z_min_layer) const { int z_coord = Kokkos::round( (raw_temperature_data(i, 2) + deltax * layer_height * layer_counter - z_min_layer[layer_counter]) / deltax); return z_coord; } // Read data from storage, obtain melting time - double getTempCoordTM(const int i) { + double getTempCoordTM(const int i) const { double t_melting = raw_temperature_data(i, 3); return t_melting; } // Read data from storage, obtain liquidus time - double getTempCoordTL(const int i) { + double getTempCoordTL(const int i) const { double t_liquidus = raw_temperature_data(i, 4); return t_liquidus; } // Read data from storage, obtain cooling rate - double getTempCoordCR(const int i) { + double getTempCoordCR(const int i) const { double cooling_rate = raw_temperature_data(i, 5); return cooling_rate; } - // Initialize temperature fields for layer "layernumber" in case where temperature data comes from file(s) - // TODO: This can be performed on the device as the dirS problem is + // Initialize temperature fields for layer "layernumber" in case where temperature data comes from file(s), could be + // done on device but currently relies on vector/host functions as is the case for the placeNuclei function void initialize(const int layernumber, const int id, const Grid &grid, const double freezing_range, - const double deltat, const std::string simulation_type) { - + const double deltat) { + + // Counter starts at 0 for each layer + liquidus_time_counter = 0; + // Resize for new layer's domain size + Kokkos::realloc(last_time_below_liquidus, grid.domain_size); + Kokkos::realloc(current_cooling_rate, grid.domain_size); + // Init values for views updated as heating above/cooling below the liquidus occur + Kokkos::deep_copy(last_time_below_liquidus, std::numeric_limits::max()); + Kokkos::deep_copy(current_cooling_rate, 0.0); + // LayerID for the next layer + layer_id = getLayerSubview(layer_id_all_layers, grid.layer_range); // Data was already read into the "raw_temperature_data" data structure // Determine which section of "raw_temperature_data" is relevant for this layer of the overall domain int start_range = first_value[layernumber]; int end_range = last_value[layernumber]; + num_liquidus_times_this_layer = 2 * (end_range - start_range); - // Temporary host view for the maximum number of times a cell in a given layer will solidify (different for each - // layer) - view_type_int_host max_solidification_events_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), max_solidification_events); - calcMaxSolidificationEvents(id, layernumber, max_solidification_events_host, start_range, end_range, grid, - simulation_type); - int max_num_solidification_events = max_solidification_events_host(layernumber); - - // Resize liquidus_time now that the max number of solidification events is known for this layer - Kokkos::resize(liquidus_time, grid.domain_size, max_num_solidification_events, 2); - - // These views are initialized to zeros on the host, filled with data, and then copied to the device for layer - // "layernumber" - view_type_int_3d_host liquidus_time_host("LiquidusTime_H", grid.domain_size, max_num_solidification_events, 2); - view_type_float_2d_host cooling_rate_host("CoolingRate_H", grid.domain_size, max_num_solidification_events); - view_type_int_host number_of_solidification_events_host("NumSEvents_H", grid.domain_size); + // Get the number of times each cell goes below the liquidus, and the max number of solidification events in a + // cell + Kokkos::realloc(number_of_solidification_events, grid.domain_size); + number_of_solidification_events = + calcNumberOfSolidificationEvents(id, layernumber, start_range, end_range, grid); + std::vector liquidus_time_step_read(num_liquidus_times_this_layer), + cell_read(num_liquidus_times_this_layer); + std::vector cooling_rate_read(num_liquidus_times_this_layer); double largest_time = 0; double largest_time_global = 0; @@ -696,17 +769,20 @@ struct Temperature { // 1D cell coordinate on this MPI rank's domain int index = grid.get1DIndex(coord_x, coord_y, coord_z); - // Store TM, TL, CR values for this solidification event in liquidus_time_host/cooling_rate_host - const int n_solidification_events_cell = number_of_solidification_events_host(index); - liquidus_time_host(index, n_solidification_events_cell, 0) = Kokkos::round(t_melting / deltat) + 1; - liquidus_time_host(index, n_solidification_events_cell, 1) = Kokkos::round(t_liquidus / deltat) + 1; - // Cannot go above/below liquidus on same time step - must offset by at least 1 - if (liquidus_time_host(index, n_solidification_events_cell, 0) == - liquidus_time_host(index, n_solidification_events_cell, 1)) - liquidus_time_host(index, n_solidification_events_cell, 1)++; - cooling_rate_host(index, n_solidification_events_cell) = std::abs(cooling_rate_cell) * deltat; - // Increment number of solidification events for this cell - number_of_solidification_events_host(index)++; + + // Two events for this cell - going above and then going below liquidus temperature + int idx_this_layer = 2 * (i - start_range); + liquidus_time_step_read[idx_this_layer] = Kokkos::round(t_melting / deltat) + 1; + cell_read[idx_this_layer] = index; + cooling_rate_read[idx_this_layer] = -1.0; // heating + + liquidus_time_step_read[idx_this_layer + 1] = Kokkos::round(t_liquidus / deltat) + 1; + // Ensure that cell doesn't go above and below liquidus on the same time step + if (liquidus_time_step_read[idx_this_layer + 1] == liquidus_time_step_read[idx_this_layer]) + liquidus_time_step_read[idx_this_layer + 1] = liquidus_time_step_read[idx_this_layer + 1] + 1; + cell_read[idx_this_layer + 1] = index; + cooling_rate_read[idx_this_layer + 1] = std::abs(cooling_rate_cell) * deltat; // cooling + // Estimate of the time step where the last possible solidification is expected to occur double solidus_time = t_liquidus + freezing_range / cooling_rate_cell; if (solidus_time > largest_time) @@ -721,66 +797,21 @@ struct Temperature { if (id == 0) std::cout << "Layer " << layernumber << " temperatures read" << std::endl; - // Reorder solidification events in liquidus_time_host(location,event number,component) and - // cooling_rate_host(location,event number) so that they are in order based on the melting time values - // (component = 0 in liquidus_time_host) - for (int index = 0; index < grid.domain_size; index++) { - int n_solidification_events_cell = number_of_solidification_events_host(index); - if (n_solidification_events_cell > 0) { - for (int i = 0; i < n_solidification_events_cell - 1; i++) { - for (int j = (i + 1); j < n_solidification_events_cell; j++) { - if (liquidus_time_host(index, i, 0) > liquidus_time_host(index, j, 0)) { - // Swap these two points - melting event "j" happens before event "i" - float old_melt_val = liquidus_time_host(index, i, 0); - float old_liq_val = liquidus_time_host(index, i, 1); - float old_cr_val = cooling_rate_host(index, i); - liquidus_time_host(index, i, 0) = liquidus_time_host(index, j, 0); - liquidus_time_host(index, i, 1) = liquidus_time_host(index, j, 1); - cooling_rate_host(index, i) = cooling_rate_host(index, j); - liquidus_time_host(index, j, 0) = old_melt_val; - liquidus_time_host(index, j, 1) = old_liq_val; - cooling_rate_host(index, j) = old_cr_val; - } - } - } - } - } - // If a cell melts twice before reaching the liquidus temperature, this is a double counted solidification - // event and should be removed - for (int index = 0; index < grid.domain_size; index++) { - int n_solidification_events_cell = number_of_solidification_events_host(index); - if (n_solidification_events_cell > 1) { - for (int i = 0; i < n_solidification_events_cell - 1; i++) { - if (liquidus_time_host(index, i + 1, 0) < liquidus_time_host(index, i, 1)) { - std::cout << "Cell " << index << " removing anomalous event " << i + 1 << " out of " - << n_solidification_events_cell - 1 << std::endl; - // Keep whichever event has the larger liquidus time - if (liquidus_time_host(index, i + 1, 1) > liquidus_time_host(index, i, 1)) { - liquidus_time_host(index, i, 0) = liquidus_time_host(index, i + 1, 0); - liquidus_time_host(index, i, 1) = liquidus_time_host(index, i + 1, 1); - cooling_rate_host(index, i) = cooling_rate_host(index, i + 1); - } - liquidus_time_host(index, i + 1, 0) = 0.0; - liquidus_time_host(index, i + 1, 1) = 0.0; - cooling_rate_host(index, i + 1) = 0.0; - // Reshuffle other solidification events over if needed - for (int ii = (i + 1); ii < n_solidification_events_cell - 1; ii++) { - liquidus_time_host(index, ii, 0) = liquidus_time_host(index, ii + 1, 0); - liquidus_time_host(index, ii, 1) = liquidus_time_host(index, ii + 1, 1); - cooling_rate_host(index, ii) = cooling_rate_host(index, ii + 1); - } - number_of_solidification_events_host(index)--; - } - } - } - } + // Reorder liquidus time events based on the times at which they occur + initOrderedTimeTempHistory(liquidus_time_step_read, cell_read, cooling_rate_read, + num_liquidus_times_this_layer); - // Copy host view data back to device - max_solidification_events = Kokkos::create_mirror_view_and_copy(memory_space(), max_solidification_events_host); - liquidus_time = Kokkos::create_mirror_view_and_copy(memory_space(), liquidus_time_host); - cooling_rate = Kokkos::create_mirror_view_and_copy(memory_space(), cooling_rate_host); - number_of_solidification_events = - Kokkos::create_mirror_view_and_copy(memory_space(), number_of_solidification_events_host); + // LayerID data on device for this layer - assign the current LayerID if this cell appears at least once on the + // list of melt/solidification events + auto _layer_id = layer_id; + auto _liquidus_cell_list = liquidus_cell_list; + auto policy = Kokkos::RangePolicy(0, num_liquidus_times_this_layer); + Kokkos::parallel_for( + "LayerIDInit", policy, KOKKOS_LAMBDA(const int &event_num) { + const int index = _liquidus_cell_list(event_num); + _layer_id(index) = layernumber; + }); + Kokkos::fence(); MPI_Barrier(MPI_COMM_WORLD); if (id == 0) { std::cout << "Layer " << layernumber << " temperature field is from Z = " << grid.z_layer_bottom @@ -789,19 +820,10 @@ struct Temperature { } } - // Get the subview associated with the undercooling of cells in the current layer. Do not reset the undercooling of - // cells from the prior layer to zero as this information will be stored for a potential print (and a cell that - // remelts in the current layer will have its undercooling reset to 0 and recalculated) - void getCurrentLayerUndercooling(std::pair layer_range) { - undercooling_current = Kokkos::subview(undercooling_current_all_layers, layer_range); - } - - // (Optional based on selected inputs) Get the subview associated with the initial undercooling of cells during - // solidification start in the current layer. Do not reset the undercooling of cells from the prior layer to zero as - // this information will be stored for a potential print (and a cell that remelts in the current layer will have its - // undercooling reset to 0 and recalculated) - void getCurrentLayerStartingUndercooling(std::pair layer_range) { - undercooling_solidification_start = Kokkos::subview(undercooling_solidification_start_all_layers, layer_range); + // Get the subview associated with data from one layer of a multilayer simulation + template + return_view_type getLayerSubview(return_view_type input_view, std::pair layer_range) { + return Kokkos::subview(input_view, layer_range); } // For each Z coordinate, find the smallest undercooling at which solidification started and finished, writing this @@ -858,144 +880,104 @@ struct Temperature { return start_end_solidification_z_host; } - // Reset local cell undercooling (and if needed, the cell's starting undercooling) to 0 - KOKKOS_INLINE_FUNCTION - void resetUndercooling(const int index) const { - if (_store_solidification_start) - undercooling_solidification_start(index) = 0.0; - undercooling_current(index) = 0.0; - } - - // Update local cell undercooling for the current melt-resolidification event + // (Optional based on inputs) Set the starting undercooling in the cell for the solidification event that just + // started. If the cell is still above the liquid, set value to 0 KOKKOS_INLINE_FUNCTION - void updateUndercooling(const int index) const { - undercooling_current(index) += cooling_rate(index, solidification_event_counter(index)); + void setStartingUndercooling(const int cycle, const int index) const { + if (store_undercooling_solidification_start) + undercooling_solidification_start(index) = Kokkos::fmin( + 0.0, static_cast(cycle - last_time_below_liquidus(index)) * current_cooling_rate(index)); } - // (Optional based on inputs) Set the starting undercooling in the cell for the solidification event that just - // started + // (Optional based on inputs) Set the undercooling in the cell at which solidification completed KOKKOS_INLINE_FUNCTION - void setStartingUndercooling(const int index) const { - if (_store_solidification_start) - undercooling_solidification_start(index) = undercooling_current(index); + void setEndingUndercooling(const int cycle, const int index) const { + if (store_undercooling_current) + undercooling_current(index) = + static_cast(cycle - last_time_below_liquidus(index)) * current_cooling_rate(index); } // Update the solidification event counter for a cell that has not finished the previous solidification event (i.e., // solidification is not complete in this cell as it is not tempsolid or solid type) KOKKOS_INLINE_FUNCTION - void updateSolidificationCounter(const int index) const { solidification_event_counter(index)++; } - - // Update the solidification event counter for the cell (which is either tempsolid or solid type) and return whether - // all solidification events have completed in the cell - KOKKOS_INLINE_FUNCTION - bool updateCheckSolidificationCounter(const int index) const { - bool solidification_complete; - solidification_event_counter(index)++; - if (solidification_event_counter(index) == number_of_solidification_events(index)) - solidification_complete = true; - else - solidification_complete = false; - return solidification_complete; + void updateSolidificationCounter(const int index) const { + if (store_solidification_event_counter) + solidification_event_counter(index)++; } - // Reset solidification event counter and get the subview associated with the undercooling field for the next layer + // For the optionally stored views, reset solidification event counter and get the subview associated with the + // undercooling field for the next layer void resetLayerEventsUndercooling(const Grid &grid) { - getCurrentLayerUndercooling(grid.layer_range); - if (_store_solidification_start) - getCurrentLayerStartingUndercooling(grid.layer_range); - Kokkos::realloc(solidification_event_counter, grid.domain_size); - Kokkos::deep_copy(solidification_event_counter, 0); - } - - // Extract the next time that this point undergoes melting - KOKKOS_INLINE_FUNCTION - int getMeltTimeStep(const int cycle, const int index) const { - int melt_time_step; - int solidification_event_counter_cell = solidification_event_counter(index); - melt_time_step = liquidus_time(index, solidification_event_counter_cell, 0); - if (cycle > melt_time_step) { - // If the cell has already exceeded the melt time step for the current melt-solidification event, get the - // melt time step associated with the next solidification event - or, if there is no next - // melt-solidification event, return the max possible int as the cell will not melt again during this layer - // of the multilayer problem - if (solidification_event_counter_cell < (number_of_solidification_events(index) - 1)) - melt_time_step = liquidus_time(index, solidification_event_counter_cell + 1, 0); - else - melt_time_step = INT_MAX; + if (store_undercooling_current) + undercooling_current = getLayerSubview(undercooling_current_all_layers, grid.layer_range); + if (store_undercooling_solidification_start) + undercooling_solidification_start = + getLayerSubview(undercooling_solidification_start_all_layers, grid.layer_range); + if (store_solidification_event_counter) { + Kokkos::realloc(solidification_event_counter, grid.domain_size); + Kokkos::deep_copy(solidification_event_counter, 0); } - return melt_time_step; } - // Extract the next time that this point cools below the liquidus - // Uses the current value of the solidification event counter - KOKKOS_INLINE_FUNCTION - int getCritTimeStep(const int index) const { - int solidification_event_counter_cell = solidification_event_counter(index); - int crit_time_step = liquidus_time(index, solidification_event_counter_cell, 1); - return crit_time_step; - } - // Uses a specified solidification event + // Get undercooling of an active cell (init_undercooling is undercooling of domains with a uniform undercooling and + // no cooling rate) KOKKOS_INLINE_FUNCTION - int getCritTimeStep(const int index, const int solidification_event_counter_cell) const { - int crit_time_step = liquidus_time(index, solidification_event_counter_cell, 1); - return crit_time_step; + float getUndercooling(const int cycle, const int index) const { + return init_undercooling + + static_cast(cycle - last_time_below_liquidus(index)) * current_cooling_rate(index); } - // Extract the cooling rate associated with a specified solidificaiton event - KOKKOS_INLINE_FUNCTION - float getUndercoolingChange(const int index, const int solidification_event_counter_cell) const { - float undercooling_change = cooling_rate(index, solidification_event_counter_cell); - return undercooling_change; + // Extract the last time step that each cell undergoes melting in the layer. If the cell does not undergo + // solidification, either print -1 or the specified default value + template + extracted_view_data_type extractTmData(const int domain_size, const int default_val = -1) { + extracted_view_data_type extracted_data(Kokkos::ViewAllocateWithoutInitializing("extracted_data"), domain_size); + Kokkos::deep_copy(extracted_data, default_val); + auto liquidus_cell_list_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), liquidus_cell_list); + auto cooling_rate_list_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), cooling_rate_list); + + for (int n = 0; n < num_liquidus_times_this_layer; n++) { + if (cooling_rate_list_host(n) < 0) { + const int index = liquidus_cell_list_host(n); + extracted_data(index) = liquidus_time_list_host(n); + } + } + return extracted_data; } - // Extract either the last time step that all points undergo melting in the layer or the last time they cools below - // the liquidus from liquidus_time (corresponds to solidification event number_of_solidification_events-1 - // for the cell) (can't just use subview here since number_of_solidification_events is different for each cell) If - // the cell does not undergo solidification, either print -1 or the specified default value + // Extract the last time step that each cell cools below the liquidus in the layer. If the cell does not undergo + // solidification, either print -1 or the specified default value template - extracted_view_data_type extractTmTlData(const int extracted_val, const int domain_size, - const int default_val = -1) { + extracted_view_data_type extractTlData(const int domain_size, const int default_val = -1) { extracted_view_data_type extracted_data(Kokkos::ViewAllocateWithoutInitializing("extracted_data"), domain_size); - - // Local copies for lambda capture. - auto liquidus_time_local = liquidus_time; - auto number_of_solidification_events_local = number_of_solidification_events; - - auto policy = Kokkos::RangePolicy(0, domain_size); - Kokkos::parallel_for( - "Extract_tm_tl_data", policy, KOKKOS_LAMBDA(const int &index) { - // If this cell doesn't undergo solidification at all, print -1 - const int num_solidification_events_this_cell = number_of_solidification_events_local(index); - if (num_solidification_events_this_cell == 0) - extracted_data(index) = default_val; - else - extracted_data(index) = - liquidus_time_local(index, num_solidification_events_this_cell - 1, extracted_val); - }); + Kokkos::deep_copy(extracted_data, default_val); + auto liquidus_cell_list_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), liquidus_cell_list); + auto cooling_rate_list_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), cooling_rate_list); + + for (int n = 0; n < num_liquidus_times_this_layer; n++) { + if (cooling_rate_list_host(n) >= 0) { + const int index = liquidus_cell_list_host(n); + extracted_data(index) = liquidus_time_list_host(n); + } + } return extracted_data; } - // Extract the cooling rate corresponding to solidification event number `number_of_solidification_events-1` for the - // cell (can't just use subview here since number_of_solidification_events is different for each cell) If the cell - // does not undergo solidification, either print -1 or the specified default value + // Extract the cooling rate associated with the final time the cell cools below the liquidus in the layer. If + // the cell does not undergo solidification, either print -1 or the specified default value template extracted_view_data_type extractCrData(const int domain_size, const int default_val = -1) { extracted_view_data_type extracted_data(Kokkos::ViewAllocateWithoutInitializing("extracted_data"), domain_size); - - // Local copies for lambda capture. - auto cooling_rate_local = cooling_rate; - auto number_of_solidification_events_local = number_of_solidification_events; - - auto policy = Kokkos::RangePolicy(0, domain_size); - Kokkos::parallel_for( - "Extract_tm_tl_data", policy, KOKKOS_LAMBDA(const int &index) { - // If this cell doesn't undergo solidification at all, print -1 - const int num_solidification_events_this_cell = number_of_solidification_events_local(index); - if (num_solidification_events_this_cell == 0) - extracted_data(index) = default_val; - else - extracted_data(index) = cooling_rate_local(index, num_solidification_events_this_cell - 1); - }); + Kokkos::deep_copy(extracted_data, default_val); + auto liquidus_cell_list_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), liquidus_cell_list); + auto cooling_rate_list_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), cooling_rate_list); + + for (int n = 0; n < num_liquidus_times_this_layer; n++) { + if (cooling_rate_list_host(n) >= 0) { + const int index = liquidus_cell_list_host(n); + extracted_data(index) = cooling_rate_list_host(n); + } + } return extracted_data; } }; diff --git a/src/CAtimers.hpp b/src/CAtimers.hpp index 8aebcafd..a3fd2886 100644 --- a/src/CAtimers.hpp +++ b/src/CAtimers.hpp @@ -51,7 +51,7 @@ struct Timers { int id; Timer init, run, output; - Timer nucl, create_sv, capture, ghost; + Timer nucl, melt_act, create_sv, capture, ghost; Timer layer; Timer heat_transfer; @@ -61,6 +61,7 @@ struct Timers { , run() , output() , nucl() + , melt_act() , create_sv() , capture() , ghost() @@ -85,6 +86,9 @@ struct Timers { void startNucleation() { nucl.start(); } void stopNucleation() { nucl.stop(); } + void startMeltAct() { melt_act.start(); } + void stopMeltAct() { melt_act.stop(); } + void startSV() { create_sv.start(); } void stopSV() { create_sv.stop(); } @@ -122,6 +126,8 @@ struct Timers { << "]," << std::endl; log << " \"MaxMinInitTime\": [" << init.maxTime() << "," << init.minTime() << "]," << std::endl; log << " \"MaxMinNucleationTime\": [" << nucl.maxTime() << "," << nucl.minTime() << "]," << std::endl; + log << " \"MaxMinMeltingActivationTime\": [" << melt_act.maxTime() << "," << melt_act.minTime() << "]," + << std::endl; log << " \"MaxMinSteeringVectorCreationTime\": [" << create_sv.maxTime() << "," << create_sv.minTime() << "]," << std::endl; log << " \"MaxMinCellCaptureTime\": [" << capture.maxTime() << "," << capture.minTime() << "]," @@ -137,6 +143,7 @@ struct Timers { // Reduce all times across MPI ranks init.reduceMPI(); nucl.reduceMPI(); + melt_act.reduceMPI(); create_sv.reduceMPI(); capture.reduceMPI(); ghost.reduceMPI(); @@ -160,14 +167,13 @@ struct Timers { std::cout << init.print("initializing data"); std::cout << run.print("performing CA calculations"); std::cout << output.print("collecting and printing output data"); - std::cout << init.printMinMax("initializing data"); std::cout << nucl.printMinMax("in CA nucleation"); + std::cout << melt_act.printMinMax("in CA melting and reactivation of cells"); std::cout << create_sv.printMinMax("in CA steering vector creation"); std::cout << capture.printMinMax("in CA cell capture"); std::cout << ghost.printMinMax("in CA cell communication"); std::cout << output.printMinMax("exporting data"); - std::cout << "===================================================================================" << std::endl; } }; diff --git a/src/CAupdate.hpp b/src/CAupdate.hpp index 163758e9..f9aaeea2 100644 --- a/src/CAupdate.hpp +++ b/src/CAupdate.hpp @@ -19,140 +19,176 @@ #include -// For the case where all cells solidify once, determine which cells are associated with the "steering vector" of -// cells that are either active, or becoming active this time step +// For problems where the entire domain starts as liquid/active cells and then solidifies, initialize octahedra for +// initial active cells template -void fillSteeringVector_NoRemelt(const int cycle, const Grid &grid, CellData &celldata, - Temperature &temperature, Interface &interface) { +void createOctahedra_NoRemelt(const Grid &grid, CellData &celldata, Temperature &temperature, + Orientation &orientation, Interface &interface) { + + auto grain_id = celldata.getGrainIDSubview(grid); - // Cells associated with this layer that are not solid type but have passed the liquidus (crit time step) have - // their undercooling values updated Cells that meet the aforementioned criteria and are active type should be - // added to the steering vector Kokkos::parallel_for( - "FillSV", grid.domain_size, KOKKOS_LAMBDA(const int &index) { - int cell_type = celldata.cell_type(index); - bool is_not_solid = (cell_type != Solid); - int crit_time_step = temperature.liquidus_time(index, 0, 1); - bool past_crit_time = (cycle > crit_time_step); - bool cell_active = ((cell_type == Active) || (cell_type == FutureActive)); - if (is_not_solid && past_crit_time) { - temperature.updateUndercooling(index); - if (cell_active) { - interface.steering_vector(Kokkos::atomic_fetch_add(&interface.num_steer(0), 1)) = index; - } + "InitSV", grid.domain_size, KOKKOS_LAMBDA(const int &index) { + if (celldata.cell_type(index) == Active) { + // Create octahedra for cells initially designated as active ( + const int coord_x = grid.getCoordX(index); + const int coord_y = grid.getCoordY(index); + const int coord_z = grid.getCoordZ(index); + + const int my_grain_id = grain_id(index); + // The orientation for the new grain will depend on its Grain ID + const int my_orientation = getGrainOrientation(my_grain_id, orientation.n_grain_orientations); + + temperature.setStartingUndercooling(0, index); + + // Initialize new octahedron + interface.createNewOctahedron(index, coord_x, coord_y, grid.y_offset, coord_z); + + // Octahedron center is at (cx, cy, cz) - note that the Y coordinate is relative to the domain + // origin to keep the coordinate system continuous across ranks + const float cx = coord_x + 0.5; + const float cy = coord_y + grid.y_offset + 0.5; + const float cz = coord_z + 0.5; + + interface.calcCritDiagonalLength(index, cx, cy, cz, cx, cy, cz, my_orientation, + orientation.grain_unit_vector); } }); - Kokkos::deep_copy(interface.num_steer_host, interface.num_steer); } // For the case where cells may melt and solidify multiple times, determine which cells are associated with the // "steering vector" of cells that are either active, or becoming active this time step, or undergoing melting template -void fillSteeringVector_Remelt(const int cycle, const Grid &grid, const InterfacialResponseFunction &irf, - CellData &celldata, Temperature &temperature, - Interface &interface) { +void remeltActivateCells(const int cycle, const Grid &grid, const InterfacialResponseFunction &irf, + CellData &celldata, Temperature &temperature, + Interface &interface) { auto grain_id = celldata.getGrainIDSubview(grid); auto phase_id = celldata.getPhaseIDSubview(grid); - Kokkos::parallel_for( - "FillSV_RM", grid.domain_size, KOKKOS_LAMBDA(const int &index) { - int celltype = celldata.cell_type(index); - // Only iterate over cells that are not Solid type - if (celltype != Solid) { - int melt_time_step = temperature.getMeltTimeStep(cycle, index); - int crit_time_step = temperature.getCritTimeStep(index); - bool at_melt_time = (cycle == melt_time_step); - bool at_crit_time = (cycle == crit_time_step); - bool past_crit_time = (cycle > crit_time_step); - if (at_melt_time) { - // Cell melts, undercooling is reset to 0 from the previous value, if any - celldata.cell_type(index) = Liquid; - temperature.resetUndercooling(index); - if (celltype != TempSolid) { - // This cell either hasn't started or hasn't finished the previous solidification event, but - // has undergone melting - increment the solidification counter to move on to the next - // melt-solidification event - temperature.updateSolidificationCounter(index); - } - // Any adjacent active cells should also be remelted, as these cells are more likely heating up - // than cooling down These are converted to the temporary FutureLiquid state, to be later - // iterated over and loaded into the steering vector as necessary Get the x, y, z coordinates of - // the cell on this MPI rank - int coord_x = grid.getCoordX(index); - int coord_y = grid.getCoordY(index); - int coord_z = grid.getCoordZ(index); - for (int l = 0; l < 26; l++) { - // "l" correpsponds to the specific neighboring cell - // Local coordinates of adjacent cell center - int neighbor_coord_x = coord_x + interface.neighbor_x[l]; - int neighbor_coord_y = coord_y + interface.neighbor_y[l]; - int neighbor_coord_z = coord_z + interface.neighbor_z[l]; - const int neighbor_index = - grid.getNeighbor1DIndex(neighbor_coord_x, neighbor_coord_y, neighbor_coord_z); - if (neighbor_index != -1) { - if (celldata.cell_type(neighbor_index) == Active) { - // Mark adjacent active cells to this as cells that should be converted into liquid, - // as they are more likely heating than cooling - celldata.cell_type(neighbor_index) = FutureLiquid; - interface.steering_vector(Kokkos::atomic_fetch_add(&interface.num_steer(0), 1)) = - neighbor_index; + // Do any cells go above/below the liquidus time on this time step? + if (temperature.liquidus_time_counter < temperature.num_liquidus_times_this_layer) { + if (cycle == temperature.liquidus_time_list_host(temperature.liquidus_time_counter)) { + bool melt_sol_check = true; + // At least one cell goes above/below the liquidus this time step + const int first_event = temperature.liquidus_time_counter; + // Are there any other events this time step to check? + while (melt_sol_check) { + temperature.liquidus_time_counter++; + // If the previous nucleation event was the last one for this layer of the simulation, exit loop + if (temperature.liquidus_time_counter == temperature.num_liquidus_times_this_layer) + break; + // If the next nucleation event corresponds to a future time step, finish check + if (cycle != temperature.liquidus_time_list_host(temperature.liquidus_time_counter)) + melt_sol_check = false; + } + const int last_event = temperature.liquidus_time_counter; + // Loop over each event, perform associated cell type transitions and steering vector additions + auto policy = Kokkos::RangePolicy<>(first_event, last_event); + Kokkos::parallel_for( + "LiquidusTimeCheck", policy, KOKKOS_LAMBDA(const int liquidus_time_counter) { + const int index = temperature.liquidus_cell_list(liquidus_time_counter); + const int celltype = celldata.cell_type(index); + const bool cooling_yn = (temperature.cooling_rate_list(liquidus_time_counter) >= 0); + if (cooling_yn) { + temperature.last_time_below_liquidus(index) = cycle; + temperature.current_cooling_rate(index) = temperature.cooling_rate_list(liquidus_time_counter); + if ((celltype == Liquid) && (grain_id(index) != 0)) { + // Get the x, y, z coordinates of the cell on this MPI rank + int coord_x = grid.getCoordX(index); + int coord_y = grid.getCoordY(index); + int coord_z = grid.getCoordZ(index); + // If this cell has cooled to the liquidus temperature, borders at least one solid + // cell, and is part of a grain, it should become active. This only needs to be checked on + // the time step where the cell reaches the liquidus, not every time step beyond this + for (int l = 0; l < 26; l++) { + // "l" correpsponds to the specific neighboring cell + // Local coordinates of adjacent cell center + int neighbor_coord_x = coord_x + interface.neighbor_x[l]; + int neighbor_coord_y = coord_y + interface.neighbor_y[l]; + int neighbor_coord_z = coord_z + interface.neighbor_z[l]; + const int neighbor_index = + grid.getNeighbor1DIndex(neighbor_coord_x, neighbor_coord_y, neighbor_coord_z); + if (neighbor_index != -1) { + // TODO: Should check if at global domain edge in X or Y too, since this would also + // make the cell a candidate for activation. This check could also be done outside + // of the loop over neighbors l=0:25 since coord_z does not vary inside this loop + if ((celldata.cell_type(neighbor_index) == Solid) || (coord_z == 0)) { + // Cell activation to be performed as part of steering vector + l = 26; + interface.steering_vector( + Kokkos::atomic_fetch_add(&interface.num_steer(0), 1)) = index; + celldata.cell_type(index) = FutureActive; + // Single phase alloy or second phase transformation after the first + // solidification event, set phase to primary. Otherwise, set phase to secondary + const int phase_id_old = phase_id(index); + phase_id(index) = irf.getPreferredPhaseActivation(phase_id_old); + // This cell was at the edge of the temperature field - set indicator to true if + // this is being tracked + celldata.setMeltEdge(index, true); + } + } } } } - } - else if ((celltype != TempSolid) && (past_crit_time)) { - // Update cell undercooling - temperature.updateUndercooling(index); - if (celltype == Active) { - // Add active cells below liquidus to steering vector - interface.steering_vector(Kokkos::atomic_fetch_add(&interface.num_steer(0), 1)) = index; - } - } - else if ((at_crit_time) && (celltype == Liquid) && (grain_id(index) != 0)) { - // Get the x, y, z coordinates of the cell on this MPI rank - int coord_x = grid.getCoordX(index); - int coord_y = grid.getCoordY(index); - int coord_z = grid.getCoordZ(index); - // If this cell has cooled to the liquidus temperature, borders at least one solid/tempsolid - // cell, and is part of a grain, it should become active. This only needs to be checked on the - // time step where the cell reaches the liquidus, not every time step beyond this - for (int l = 0; l < 26; l++) { - // "l" correpsponds to the specific neighboring cell - // Local coordinates of adjacent cell center - int neighbor_coord_x = coord_x + interface.neighbor_x[l]; - int neighbor_coord_y = coord_y + interface.neighbor_y[l]; - int neighbor_coord_z = coord_z + interface.neighbor_z[l]; - const int neighbor_index = - grid.getNeighbor1DIndex(neighbor_coord_x, neighbor_coord_y, neighbor_coord_z); - if (neighbor_index != -1) { - if ((celldata.cell_type(neighbor_index) == TempSolid) || - (celldata.cell_type(neighbor_index) == Solid) || (coord_z == 0)) { - // Cell activation to be performed as part of steering vector - l = 26; - interface.steering_vector(Kokkos::atomic_fetch_add(&interface.num_steer(0), 1)) = index; - celldata.cell_type(index) = FutureActive; - // Single phase alloy or second phase transformation after the first solidification - // event, set phase to primary. Otherwise, set phase to secondary - const int phase_id_old = phase_id(index); - phase_id(index) = irf.getPreferredPhaseActivation(phase_id_old); - // This cell was at the edge of the temperature field - set indicator to true if this is - // being tracked - celldata.setMeltEdge(index, true); + else { + // Cell melts, undercooling is reset to 0 from the previous value, if any + celldata.cell_type(index) = Liquid; + // Reset this to max value as this cell is no longer below the liquidus + temperature.last_time_below_liquidus(index) = std::numeric_limits::max(); + // Any adjacent active cells should also be remelted, as these cells are more likely heating up + // than cooling down These are converted to the temporary FutureLiquid state, to be later + // iterated over and loaded into the steering vector as necessary Get the x, y, z coordinates of + // the cell on this MPI rank + int coord_x = grid.getCoordX(index); + int coord_y = grid.getCoordY(index); + int coord_z = grid.getCoordZ(index); + for (int l = 0; l < 26; l++) { + // "l" correpsponds to the specific neighboring cell + // Local coordinates of adjacent cell center + int neighbor_coord_x = coord_x + interface.neighbor_x[l]; + int neighbor_coord_y = coord_y + interface.neighbor_y[l]; + int neighbor_coord_z = coord_z + interface.neighbor_z[l]; + const int neighbor_index = + grid.getNeighbor1DIndex(neighbor_coord_x, neighbor_coord_y, neighbor_coord_z); + if (neighbor_index != -1) { + if (celldata.cell_type(neighbor_index) == Active) { + // Mark adjacent active cells to this as cells that should be converted into liquid, + // as they are more likely heating than cooling + celldata.cell_type(neighbor_index) = FutureLiquid; + interface.steering_vector(Kokkos::atomic_fetch_add(&interface.num_steer(0), 1)) = + neighbor_index; + } } } } - } + }); + } + } + Kokkos::fence(); +} + +// Determine which cells are associated with the "steering vector" of cells that are active and below the liquidus +// temperature on this time step +template +void fillSteeringVector(const int cycle, const Grid &grid, CellData &celldata, + Temperature &temperature, Interface &interface) { + Kokkos::parallel_for( + "FillSV", grid.domain_size, KOKKOS_LAMBDA(const int &index) { + if ((celldata.cell_type(index) == Active) && (cycle > temperature.last_time_below_liquidus(index))) { + // Add active cells below liquidus to steering vector + interface.steering_vector(Kokkos::atomic_fetch_add(&interface.num_steer(0), 1)) = index; } }); Kokkos::fence(); - // Copy size of steering vector (containing positions of undercooled liquid/active cells) to the host + // Copy size of steering vector (containing positions of undercooled active cells added in this routine, as well as + // FutureLiquid and FutureActive cells added in remeltActivateCells) to the host Kokkos::deep_copy(interface.num_steer_host, interface.num_steer); } // Decentered octahedron algorithm for the capture of new interface cells by grains template -void cellCapture(const int, const int np, const Grid &grid, const InterfacialResponseFunction &irf, +void cellCapture(const int cycle, const int np, const Grid &grid, const InterfacialResponseFunction &irf, CellData &celldata, Temperature &temperature, Interface &interface, Orientation &orientation) { @@ -175,8 +211,9 @@ void cellCapture(const int, const int np, const Grid &grid, const InterfacialRes // Cells of interest for the CA - active cells and future active/liquid cells if (cell_type_old == Active) { // Get undercooling of active cell - const float local_undercooling = temperature.undercooling_current(index); - // Update diagonal length of octahedron based on local undercooling and interfacial response function + const float local_undercooling = temperature.getUndercooling(cycle, index); + // Update diagonal length of octahedron based on local undercooling and interfacial response + // function interface.diagonal_length(index) += irf.compute(local_undercooling, phase_id(index)); const float diagonal_length_cell = interface.diagonal_length(index); // Switch that becomes false if the cell has at least 1 liquid type neighbor @@ -225,7 +262,7 @@ void cellCapture(const int, const int np, const Grid &grid, const InterfacialRes phase_id(neighbor_index) = my_phase_id; // Store the initial undercooling for the newly captured cell, if this output was // toggled - temperature.setStartingUndercooling(neighbor_index); + temperature.setStartingUndercooling(cycle, neighbor_index); // (cxold, cyold, czold) are the coordinates of this decentered octahedron const float cxold = interface.octahedron_center(3 * index); @@ -440,16 +477,12 @@ void cellCapture(const int, const int np, const Grid &grid, const InterfacialRes } // End if statement over neighbors on the active grid } // End loop over all neighbors of this active cell if (deactivate_cell) { - // This active cell has no more neighboring cells to be captured - // Update the counter for the number of times this cell went from liquid to active to solid - bool solidification_complete_y_n = temperature.updateCheckSolidificationCounter(index); - // Did the cell solidify for the last time in the layer? - // If so, this cell is solid - ignore until next layer (if needed) - // If not, this cell is tempsolid, will become liquid again - if (solidification_complete_y_n) - celldata.cell_type(index) = Solid; - else - celldata.cell_type(index) = TempSolid; + // This active cell has no more neighboring cells to be captured, is now solid + celldata.cell_type(index) = Solid; + // Set undercooling of cell at solidification completion/update counter if optional inputs were + // toggled + temperature.setEndingUndercooling(cycle, index); + temperature.updateSolidificationCounter(index); } } else if (cell_type_old == FutureActive) { @@ -459,6 +492,8 @@ void cellCapture(const int, const int np, const Grid &grid, const InterfacialRes const int my_grain_id = grain_id(index); // grain_id was assigned as part of Nucleation const int my_phase_id = phase_id(index); // phase_id was assigned as part of Nucleation + temperature.setStartingUndercooling(cycle, index); + // Initialize new octahedron interface.createNewOctahedron(index, coord_x, coord_y, grid.y_offset, coord_z); // The orientation for the new grain will depend on its Grain ID (nucleated grains have negative @@ -743,52 +778,33 @@ void haloUpdate(const int, const int, const Grid &grid, CellData &c // Jump to the next time step with work to be done, if no melting or solidification events occur in the next 5000 time // steps template -void jumpTimeStep(int &cycle, int remaining_liquid_cells, const int local_temp_solid_cells, - Nucleation &nucleation, Temperature &temperature, const Grid &grid, - CellData &celldata, const int id, const int layernumber, const int np, +void jumpTimeStep(int &cycle, Nucleation &nucleation, Temperature &temperature, + const Grid &grid, CellData &celldata, const int id, const int layernumber, const int np, Orientation &orientation, Print &print, const double deltat, Interface &interface) { - MPI_Bcast(&remaining_liquid_cells, 1, MPI_INT, 0, MPI_COMM_WORLD); - if (remaining_liquid_cells == 0) { - // If this rank still has cells that will later undergo transformation (TempSolid), check when - // the next solid cells go above the liquidus (remelting) Otherwise, assign the largest possible time step as - // the next time work needs to be done on the rank - int next_melt_time_step; - if (local_temp_solid_cells > 0) { - Kokkos::parallel_reduce( - "CheckNextTSForWork", grid.domain_size, - KOKKOS_LAMBDA(const int &index, int &tempv) { - if (celldata.cell_type(index) == TempSolid) { - int solidification_counter_this_cell = temperature.solidification_event_counter(index); - int next_melt_time_step_this_cell = - temperature.liquidus_time(index, solidification_counter_this_cell, 0); - if (next_melt_time_step_this_cell < tempv) - tempv = next_melt_time_step_this_cell; - } - }, - Kokkos::Min(next_melt_time_step)); - } - else - next_melt_time_step = INT_MAX; - - int global_next_melt_time_step; - MPI_Allreduce(&next_melt_time_step, &global_next_melt_time_step, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); - if ((global_next_melt_time_step - cycle) > 5000) { - // Print current state of the system for desired output fields for any of the time - // steps between now and when melting/solidification occurs again, if the print option for idle frame - // printing was toggled - - print.printIdleIntralayer(id, np, layernumber, deltat, cycle, grid, celldata, temperature, interface, - orientation, global_next_melt_time_step); - // Jump to next time step when melting occurs - cycle = global_next_melt_time_step - 1; - // If nucleation events were possible in any of the skipped time steps, update the counter accordingly as - // these nucleation events could not have occured without any available liquid cells - nucleation.advanceCounterSkippedTimeSteps(cycle); - if (id == 0) - std::cout << "Jumping to cycle " << cycle + 1 << std::endl; - } + int local_next_melt_time_step; + if (temperature.liquidus_time_counter == temperature.num_liquidus_times_this_layer) + local_next_melt_time_step = std::numeric_limits::max(); + else + local_next_melt_time_step = temperature.liquidus_time_list_host(temperature.liquidus_time_counter); + int next_melt_time_step; + MPI_Allreduce(&local_next_melt_time_step, &next_melt_time_step, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); + + if ((next_melt_time_step - cycle) > 5000) { + // Print current state of the system for desired output fields for any of the time + // steps between now and when melting/solidification occurs again, if the print option for idle frame + // printing was toggled + + print.printIdleIntralayer(id, np, layernumber, deltat, cycle, grid, celldata, temperature, interface, + orientation, next_melt_time_step); + // Jump to next time step when melting occurs + cycle = next_melt_time_step - 1; + // If nucleation events were possible in any of the skipped time steps, update the counter accordingly as + // these nucleation events could not have occured without any available liquid cells + nucleation.advanceCounterSkippedTimeSteps(cycle); + if (id == 0) + std::cout << "Jumping to cycle " << cycle + 1 << std::endl; } } @@ -803,40 +819,34 @@ void intermediateOutputAndCheck(const int id, const int np, int &cycle, const Gr Interface &interface) { auto grain_id = celldata.getGrainIDSubview(grid); - int local_superheated_cells, local_undercooled_cells, local_active_cells, local_temp_solid_cells, - local_finished_solid_cells; + int local_superheated_cells, local_undercooled_cells, local_active_cells, local_solid_cells; Kokkos::parallel_reduce( "IntermediateOutput", grid.domain_size, - KOKKOS_LAMBDA(const int &index, int &sum_superheated, int &sum_undercooled, int &sum_active, - int &sum_temp_solid, int &sum_finished_solid) { + KOKKOS_LAMBDA(const int &index, int &sum_superheated, int &sum_undercooled, int &sum_active, int &sum_solid) { int cell_type_this_cell = celldata.cell_type(index); if (cell_type_this_cell == Liquid) { - int crit_time_step = temperature.getCritTimeStep(index); - if (crit_time_step > cycle) + if (cycle < temperature.last_time_below_liquidus(index)) sum_superheated += 1; else sum_undercooled += 1; } else if (cell_type_this_cell == Active) sum_active += 1; - else if (cell_type_this_cell == TempSolid) - sum_temp_solid += 1; else if (cell_type_this_cell == Solid) - sum_finished_solid += 1; + sum_solid += 1; }, - local_superheated_cells, local_undercooled_cells, local_active_cells, local_temp_solid_cells, - local_finished_solid_cells); + local_superheated_cells, local_undercooled_cells, local_active_cells, local_solid_cells); int global_successful_nuc_events_this_rank = 0; - int global_superheated_cells, global_undercooled_cells, global_active_cells, global_temp_solid_cells, - global_finished_solid_cells; + int global_superheated_cells, global_undercooled_cells, global_active_cells, global_solid_cells; MPI_Reduce(&local_superheated_cells, &global_superheated_cells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(&local_undercooled_cells, &global_undercooled_cells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(&local_active_cells, &global_active_cells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); - MPI_Reduce(&local_temp_solid_cells, &global_temp_solid_cells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); - MPI_Reduce(&local_finished_solid_cells, &global_finished_solid_cells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + MPI_Reduce(&local_solid_cells, &global_solid_cells, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); MPI_Reduce(&successful_nuc_events_this_rank, &global_successful_nuc_events_this_rank, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); + // Cells of interest are those currently undergoing a melting-solidification cycle + int remaining_cells_of_interest; if (id == 0) { std::cout << "Current time step " << cycle << " on layer number " << layernumber << std::endl; @@ -844,21 +854,26 @@ void intermediateOutputAndCheck(const int id, const int np, int &cycle, const Gr << "/" << global_undercooled_cells << std::endl; std::cout << "Number of active (solid-liquid interface) cells in simulation: " << global_active_cells << std::endl; - std::cout << "Number of solid cells in simulation (finished/to be remelted): " << global_finished_solid_cells - << "/" << global_temp_solid_cells << std::endl; + std::cout << "Number of solid cells in simulation (finished/to be remelted): " << global_solid_cells + << std::endl; std::cout << "Number of nucleation events during simulation of this layer: " << global_successful_nuc_events_this_rank << std::endl; std::cout << "======================================================================================" << std::endl; - if (global_superheated_cells + global_undercooled_cells + global_active_cells + global_temp_solid_cells == 0) - x_switch = 1; + remaining_cells_of_interest = global_active_cells + global_undercooled_cells + global_superheated_cells; } - MPI_Bcast(&x_switch, 1, MPI_INT, 0, MPI_COMM_WORLD); - // Cells of interest are those currently undergoing a melting-solidification cycle - int remaining_cells_of_interest = global_active_cells + global_superheated_cells + global_undercooled_cells; - if ((x_switch == 0) && ((simulation_type != "Directional"))) - jumpTimeStep(cycle, remaining_cells_of_interest, local_temp_solid_cells, nucleation, temperature, grid, - celldata, id, layernumber, np, orientation, print, deltat, interface); + int local_liquidus_events_remaining = temperature.num_liquidus_times_this_layer - temperature.liquidus_time_counter; + int liquidus_events_remaining; + MPI_Allreduce(&local_liquidus_events_remaining, &liquidus_events_remaining, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&local_liquidus_events_remaining, &liquidus_events_remaining, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + MPI_Bcast(&remaining_cells_of_interest, 1, MPI_INT, 0, MPI_COMM_WORLD); + + // Every cell has crossed the liquidus for the final time and all cells are solid + if ((liquidus_events_remaining == 0) && (remaining_cells_of_interest == 0)) + x_switch = 1; + if ((x_switch == 0) && (simulation_type != "Directional") && (remaining_cells_of_interest == 0)) + jumpTimeStep(cycle, nucleation, temperature, grid, celldata, id, layernumber, np, orientation, print, deltat, + interface); } //*****************************************************************************/ diff --git a/src/runCA.hpp b/src/runCA.hpp index 646662f4..c5afd158 100644 --- a/src/runCA.hpp +++ b/src/runCA.hpp @@ -20,20 +20,28 @@ void runExaCA(int id, int np, Inputs inputs, Timers timers, Grid grid, Temperatu using memory_space = MemorySpace; std::string simulation_type = inputs.simulation_type; + // Full domain solidification - all cells initially liquid or active and end up solid + bool full_domain_solidification; + if ((simulation_type == "Directional") || (simulation_type == "SingleGrain")) + full_domain_solidification = true; + else + full_domain_solidification = false; // Material response function InterfacialResponseFunction irf(inputs.domain.deltat, grid.deltax, inputs.irf); - // Read temperature data if necessary + // Read temperature data if necessary. For the spot problem, store spot melt data as if it were read from a file if (simulation_type == "FromFile") temperature.readTemperatureData(id, grid, 0); - // Initialize the temperature fields for the simulation type of interest - if ((simulation_type == "Directional") || (simulation_type == "SingleGrain")) - temperature.initialize(id, simulation_type, grid, inputs.domain.deltat); else if (simulation_type == "Spot") - temperature.initialize(id, grid, irf.freezingRange(), inputs.domain.deltat, inputs.domain.spot_radius); - else if ((simulation_type == "FromFile") || (simulation_type == "FromFinch")) - temperature.initialize(0, id, grid, irf.freezingRange(), inputs.domain.deltat, simulation_type); + temperature.storeSpotData(id, grid, irf.freezingRange(), inputs.domain.deltat, inputs.domain.spot_radius); + + // Initialize the temperature fields for the simulation type of interest. These are either simple unidirectional + // fields, or more complex data stored in a view in the temperature struct + if (full_domain_solidification) + temperature.initialize(id, simulation_type, grid, inputs.domain.deltat); + else + temperature.initialize(0, id, grid, irf.freezingRange(), inputs.domain.deltat); MPI_Barrier(MPI_COMM_WORLD); // Initialize grain orientations @@ -44,16 +52,19 @@ void runExaCA(int id, int np, Inputs inputs, Timers timers, Grid grid, Temperatu // Initialize cell types, grain IDs, and layer IDs CellData celldata(grid, inputs.substrate, inputs.print.store_melt_pool_edge); if (simulation_type == "Directional") - celldata.initSubstrate(id, grid, inputs.rng_seed); + celldata.initSubstrate_Directional(id, grid, inputs.rng_seed); else if (simulation_type == "SingleGrain") - celldata.initSubstrate(id, grid); + celldata.initSubstrate_SingleGrain(id, grid); else - celldata.initSubstrate(id, grid, inputs.rng_seed, temperature.number_of_solidification_events); + celldata.initSubstrate_BaseplatePowder(id, grid, inputs.rng_seed); MPI_Barrier(MPI_COMM_WORLD); // Variables characterizing the active cell region within each rank's grid, including buffers for ghost node data // (fixed size) and the steering vector/steering vector size on host/device Interface interface(id, grid.domain_size, inputs.substrate.init_oct_size); + // Initialize octahedra for initial active cells, if necessary for this problem type + if (full_domain_solidification) + createOctahedra_NoRemelt(grid, celldata, temperature, orientation, interface); MPI_Barrier(MPI_COMM_WORLD); // Nucleation data structure, containing views of nuclei locations, time steps, and ids, and nucleation event @@ -66,7 +77,7 @@ void runExaCA(int id, int np, Inputs inputs, Timers timers, Grid grid, Temperatu // Fill in nucleation data structures, and assign nucleation undercooling values to potential nucleation events // Potential nucleation grains are only associated with liquid cells in layer 0 - they will be initialized for each // successive layer when layer 0 is complete - nucleation.placeNuclei(temperature, irf, inputs.rng_seed, 0, grid, id); + nucleation.placeNuclei(simulation_type, temperature, irf, inputs.rng_seed, 0, grid, id, inputs.domain.deltat); // Initialize printing struct from inputs Print print(grid, np, inputs.print); @@ -97,14 +108,16 @@ void runExaCA(int id, int np, Inputs inputs, Timers timers, Grid grid, Temperatu nucleation.nucleateGrain(cycle, grid, celldata, interface); timers.stopNucleation(); - // Cells that have a successful nucleation event, and other cells that are at the solid-liquid interface are - // added to a steering vector. Logic in fillSteeringVector_NoRemelt is a simpified version of - // fillSteeringVector_Remelt + // Melt cells above the liquidus, activate cells at the solid-liquid interface below the liquidus + if (!full_domain_solidification) { + timers.startMeltAct(); + remeltActivateCells(cycle, grid, irf, celldata, temperature, interface); + timers.stopMeltAct(); + } + + // Create steering vector of cells that are active and undercooled on this time step timers.startSV(); - if ((simulation_type == "Directional") || (simulation_type == "SingleGrain")) - fillSteeringVector_NoRemelt(cycle, grid, celldata, temperature, interface); - else - fillSteeringVector_Remelt(cycle, grid, irf, celldata, temperature, interface); + fillSteeringVector(cycle, grid, celldata, temperature, interface); timers.stopSV(); // Iterate over the steering vector to perform active cell creation and capture operations, and if needed, @@ -156,8 +169,7 @@ void runExaCA(int id, int np, Inputs inputs, Timers timers, Grid grid, Temperatu temperature.readTemperatureData(id, grid, layernumber + 1); MPI_Barrier(MPI_COMM_WORLD); // Initialize next layer's temperature data - temperature.initialize(layernumber + 1, id, grid, irf.freezingRange(0), inputs.domain.deltat, - simulation_type); + temperature.initialize(layernumber + 1, id, grid, irf.freezingRange(0), inputs.domain.deltat); // Reset solidification event counter of all cells to zeros for the next layer, resizing to number of cells // associated with the next layer, and get the subview for undercooling @@ -168,14 +180,14 @@ void runExaCA(int id, int np, Inputs inputs, Timers timers, Grid grid, Temperatu MPI_Barrier(MPI_COMM_WORLD); // Sets up views, powder layer (if necessary), and cell types for the next layer of a multilayer problem - celldata.initNextLayer(layernumber + 1, id, grid, inputs.rng_seed, - temperature.number_of_solidification_events); + celldata.initNextLayer(layernumber + 1, id, grid, inputs.rng_seed); // Initialize potential nucleation event data for next layer "layernumber + 1" // Views containing nucleation data will be resized to the possible number of nuclei on a given MPI rank for // the next layer nucleation.resetNucleiCounters(); // start counters at 0 - nucleation.placeNuclei(temperature, irf, inputs.rng_seed, layernumber + 1, grid, id); + nucleation.placeNuclei(simulation_type, temperature, irf, inputs.rng_seed, layernumber + 1, grid, id, + inputs.domain.deltat); x_switch = 0; MPI_Barrier(MPI_COMM_WORLD); @@ -196,7 +208,7 @@ void runExaCA(int id, int np, Inputs inputs, Timers timers, Grid grid, Temperatu interface, orientation); // Calculate volume fraction of solidified domain consisting of nucleated grains - float vol_fraction_nucleated = celldata.calcVolFractionNucleated(id, grid); + float vol_fraction_nucleated = celldata.calcVolFractionNucleated(id, grid, temperature); // Get MPI timing statisticss timers.stopOutput(); timers.reduceMPI(); diff --git a/unit_test/tstCellData.hpp b/unit_test/tstCellData.hpp index 79be1f9d..873087b8 100644 --- a/unit_test/tstCellData.hpp +++ b/unit_test/tstCellData.hpp @@ -55,7 +55,7 @@ void testCellDataInit_SingleGrain() { EXPECT_EQ(inputs.substrate.single_grain_orientation, celldata._inputs.single_grain_orientation); // Init grain - celldata.initSubstrate(id, grid); + celldata.initSubstrate_SingleGrain(id, grid); // Copy cell type and grain ID back to host to check if the values match - only 1 cell should've been assigned type // active and GrainID = 1 (though it may be duplicated in the ghost nodes of other ranks) auto grain_id_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), celldata.grain_id_all_layers); @@ -68,7 +68,7 @@ void testCellDataInit_SingleGrain() { if ((coord_z == expected_grain_z) && (coord_x == expected_grain_x) && (coord_y_global == expected_grain_y)) { EXPECT_EQ(grain_id_host(coord_1d), celldata._inputs.single_grain_orientation + 1); - EXPECT_EQ(cell_type_host(coord_1d), FutureActive); + EXPECT_EQ(cell_type_host(coord_1d), Active); } else { EXPECT_EQ(grain_id_host(coord_1d), 0); @@ -122,7 +122,7 @@ void testCellDataInit_Constrained_Automatic(std::string input_surface_init_mode) // Check appropriate initialization of celldata input EXPECT_DOUBLE_EQ(inputs.substrate.fract_surface_sites_active, celldata._inputs.fract_surface_sites_active); // Initialize substrate grains - celldata.initSubstrate(id, grid, inputs.rng_seed); + celldata.initSubstrate_Directional(id, grid, inputs.rng_seed); // Copy CellType, GrainID views to host to check values auto cell_type_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), celldata.cell_type); auto grain_id_all_layers_host = @@ -142,7 +142,7 @@ void testCellDataInit_Constrained_Automatic(std::string input_surface_init_mode) else { // Check that active cells have GrainIDs > 0, and less than max_expected_grain_id (GrainIDs 1 through // max_expected_grain_id-1 should've been used) - if (cell_type_host(index) == FutureActive) { + if (cell_type_host(index) == Active) { EXPECT_GT(grain_id_all_layers_host(index), 0); EXPECT_LT(grain_id_all_layers_host(index), max_expected_grain_id); } @@ -175,13 +175,13 @@ void testCellDataInit_Constrained_Custom() { CellData celldata(grid, inputs.substrate); // Place substrate grains - celldata.initSubstrate(id, grid, 0.0); + celldata.initSubstrate_Directional(id, grid, 0.0); // Copy CellType, GrainID views to host to check values auto cell_type_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), celldata.cell_type); auto grain_id_all_layers_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), celldata.grain_id_all_layers); - // Bottom surface (Z = 0): Check that the two FutureActive cells are in the right locations on the bottom surface, + // Bottom surface (Z = 0): Check that the two Active cells are in the right locations on the bottom surface, // and have the right grain IDs. Otherwise, cells are liquid and have GrainID still equal to 0 for (int coord_x = 0; coord_x < grid.nx; coord_x++) { for (int coord_y_local = 0; coord_y_local < grid.ny_local; coord_y_local++) { @@ -189,11 +189,11 @@ void testCellDataInit_Constrained_Custom() { int coord_y_global = coord_y_local + grid.y_offset; if ((coord_x == 100) && (coord_y_global == 50)) { EXPECT_EQ(grain_id_all_layers_host(index), 25); - EXPECT_EQ(cell_type_host(index), FutureActive); + EXPECT_EQ(cell_type_host(index), Active); } else if ((coord_x == 100) && (coord_y_global == 150)) { EXPECT_EQ(grain_id_all_layers_host(index), 9936); - EXPECT_EQ(cell_type_host(index), FutureActive); + EXPECT_EQ(cell_type_host(index), Active); } else { EXPECT_EQ(grain_id_all_layers_host(index), 0); @@ -212,7 +212,6 @@ void testCellDataInit_Constrained_Custom() { void testCellDataInit(bool powder_first_layer) { using memory_space = TEST_MEMSPACE; - using view_int = Kokkos::View; int id, np; // Get number of processes @@ -275,22 +274,6 @@ void testCellDataInit(bool powder_first_layer) { inputs.substrate.baseplate_grain_spacing = 3.0; inputs.rng_seed = 0.0; - // Create dummy temperature data - view_int number_of_solidification_events(Kokkos::ViewAllocateWithoutInitializing("NumberOfSolidificationEvents"), - grid.domain_size); - Kokkos::parallel_for( - "InitTestTemperatureData", grid.domain_size, KOKKOS_LAMBDA(const int &index) { - int coord_x = grid.getCoordX(index); - int coord_y = grid.getCoordY(index); - // Assign some of these a value of 0 (these will be solid cells), and others a positive value - // (these will be tempsolid cells) - if (coord_x + coord_y % 2 == 0) - number_of_solidification_events(index) = 0; - else - number_of_solidification_events(index) = 1; - }); - Kokkos::fence(); - // Call constructor CellData celldata(grid, inputs.substrate); // Check that substrate inputs were copied from inputs struct correctly @@ -300,7 +283,7 @@ void testCellDataInit(bool powder_first_layer) { EXPECT_FALSE(celldata._inputs.baseplate_through_powder); EXPECT_DOUBLE_EQ(celldata._inputs.powder_grain_spacing, grid.deltax * pow(10, 6)); // Initialize baseplate grain structure - celldata.initSubstrate(id, grid, inputs.rng_seed, number_of_solidification_events); + celldata.initSubstrate_BaseplatePowder(id, grid, inputs.rng_seed); // Copy GrainID results back to host to check first layer's initialization auto grain_id_all_layers_host = @@ -330,22 +313,15 @@ void testCellDataInit(bool powder_first_layer) { EXPECT_EQ(celldata.next_layer_first_epitaxial_grain_id, np + 1); } - // Copy cell types back to host to check + // Copy cell types back to host to check - all cells should be initialized as solid auto cell_type_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), celldata.cell_type); - for (int index = 0; index < grid.domain_size; index++) { - int coord_x = grid.getCoordX(index); - int coord_y = grid.getCoordY(index); - // Cells with no associated solidification events should be solid, others TempSolid - if (coord_x + coord_y % 2 == 0) - EXPECT_EQ(cell_type_host(index), Solid); - else - EXPECT_EQ(cell_type_host(index), TempSolid); - } + for (int index = 0; index < grid.domain_size; index++) + EXPECT_EQ(cell_type_host(index), Solid); int previous_layer_first_epitaxial_grain_id = celldata.next_layer_first_epitaxial_grain_id; // Initialize the next layer using the same time-temperature history - powder should span cells at Z = 3 expected_num_powder_grains_per_layer = grid.nx * grid.ny_local * np; grid.initNextLayer(id, "FromFile", 1); - celldata.initNextLayer(1, id, grid, inputs.rng_seed, number_of_solidification_events); + celldata.initNextLayer(1, id, grid, inputs.rng_seed); // Copy all grain IDs for all layers back to the host to check that they match // and that the powder layer was initialized correctly @@ -372,18 +348,10 @@ void testCellDataInit(bool powder_first_layer) { EXPECT_EQ(grain_id_host(index), grain_id_all_layers_host(index_all_layers)); } - // Copy cell types back to host to check - should be the same as the previous layer as the same time-temperature - // history was used + // Copy cell types back to host to check - should have also initialized as Solid cell_type_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), celldata.cell_type); - for (int index = 0; index < grid.domain_size; index++) { - int coord_x = grid.getCoordX(index); - int coord_y = grid.getCoordY(index); - // Cells with no associated solidification events should be solid, others TempSolid - if (coord_x + coord_y % 2 == 0) - EXPECT_EQ(cell_type_host(index), Solid); - else - EXPECT_EQ(cell_type_host(index), TempSolid); - } + for (int index = 0; index < grid.domain_size; index++) + EXPECT_EQ(cell_type_host(index), Solid); } void testCalcVolFractionNucleated() { @@ -420,10 +388,8 @@ void testCalcVolFractionNucleated() { // Let all cells except those at Z = 0 have undergone solidification // Let the cells at Z = 1 consist of positive grain IDs, and those at Z = 2 of negative grain IDs CellData celldata(grid, inputs.substrate); - Kokkos::View grain_id_host(Kokkos::ViewAllocateWithoutInitializing("GrainID"), - grid.domain_size_all_layers); - Kokkos::View layer_id_host(Kokkos::ViewAllocateWithoutInitializing("LayerID"), - grid.domain_size_all_layers); + // LayerID stored in temperature view + Temperature temperature(grid, inputs.temperature, inputs.print); auto md_policy = Kokkos::MDRangePolicy>( {0, 0, 0}, {grid.nz, grid.nx, grid.ny_local}); @@ -431,9 +397,9 @@ void testCalcVolFractionNucleated() { "VolFractNucleatedInit", md_policy, KOKKOS_LAMBDA(const int coord_z, const int coord_x, const int coord_y) { int index = grid.get1DIndex(coord_x, coord_y, coord_z); if (coord_z == 0) - celldata.layer_id_all_layers(index) = -1; + temperature.layer_id_all_layers(index) = -1; else - celldata.layer_id_all_layers(index) = 0; + temperature.layer_id_all_layers(index) = 0; if (coord_z == 2) celldata.grain_id_all_layers(index) = -1; else @@ -441,7 +407,7 @@ void testCalcVolFractionNucleated() { }); // Perform calculation and compare to expected value (half of the solidified portion of the domain should consist of // nucleated grains, regardless of the number of MPI ranks used) - float vol_fraction_nucleated = celldata.calcVolFractionNucleated(id, grid); + float vol_fraction_nucleated = celldata.calcVolFractionNucleated(id, grid, temperature); EXPECT_FLOAT_EQ(vol_fraction_nucleated, 0.5); } //---------------------------------------------------------------------------// diff --git a/unit_test/tstInterface.hpp b/unit_test/tstInterface.hpp index df41b2e8..4f917096 100644 --- a/unit_test/tstInterface.hpp +++ b/unit_test/tstInterface.hpp @@ -459,11 +459,13 @@ void testResizeBuffers() { int id; MPI_Comm_rank(MPI_COMM_WORLD, &id); + Grid grid; // Interface struct int domain_size = 50; + grid.domain_size = domain_size; int buf_size_initial_estimate = 50; // Init buffers to large size - Interface interface(id, domain_size, 0.01, 50); + Interface interface(id, grid.domain_size, 0.01, 50); // Fill buffers with test data Kokkos::parallel_for( @@ -500,18 +502,21 @@ void testResizeBuffers() { //---------------------------------------------------------------------------// // Steering vector and decentered octahedron tests //---------------------------------------------------------------------------// -void testFillSteeringVector_Remelt() { +void testRemeltActivateCells() { using memory_space = TEST_MEMSPACE; + int id; + MPI_Comm_rank(MPI_COMM_WORLD, &id); + // Grid struct with manually set values - // Create views - each rank has 125 cells, 75 of which are part of the active region of the domain + // Create views - each rank has 125 cells Grid grid; grid.nx = 5; grid.ny_local = 5; grid.y_offset = 5; - grid.nz = 5; - grid.z_layer_bottom = 2; + grid.nz = 3; + grid.z_layer_bottom = 0; grid.nz_layer = 3; grid.domain_size = grid.nx * grid.ny_local * grid.nz_layer; grid.domain_size_all_layers = grid.nx * grid.ny_local * grid.nz; @@ -531,99 +536,125 @@ void testFillSteeringVector_Remelt() { inputs.irf.function[0] = inputs.irf.cubic; InterfacialResponseFunction irf(inputs.domain.deltat, grid.deltax, inputs.irf); - // Initialize cell/temperature structures + // Initialize cell/temperature/orientation structures CellData celldata(grid, inputs.substrate); auto grain_id = celldata.getGrainIDSubview(grid); - Temperature temperature(grid, inputs.temperature); - - // Fill temperature structure - Kokkos::parallel_for( - "TempDataInit", grid.domain_size, KOKKOS_LAMBDA(const int &index) { - grain_id(index) = 1; - celldata.cell_type(index) = TempSolid; - // Cell coordinates on this rank in X, Y, and Z (GlobalZ = relative to domain bottom) - int coord_z = grid.getCoordZ(index); - int coord_y = grid.getCoordY(index); - // Cells at Z = 0 through Z = 2 are Solid, Z = 3 and 4 are TempSolid - if (coord_z <= 2) { - // Solid cells should have -1 assigned as their melt/crit time steps - temperature.liquidus_time(index, 0, 0) = -1; - temperature.liquidus_time(index, 0, 1) = -1; - temperature.cooling_rate(index, 0) = 0; - } - else { - // Cells "melt" at a time step corresponding to their Y location in the overall domain (depends on - // y_offset of the rank) - temperature.liquidus_time(index, 0, 0) = coord_y + grid.y_offset + 1; - // Cells reach liquidus during cooling 2 time steps after melting - temperature.liquidus_time(index, 0, 1) = temperature.liquidus_time(index, 0, 0) + 2; - } - temperature.cooling_rate(index, 0) = 0.2; - temperature.number_of_solidification_events(index) = 1; - }); + Temperature temperature(grid, inputs.temperature, inputs.print); + std::string grain_orientation_file_tmp = checkFileInstalled("GrainOrientationVectors.csv", id); + std::vector grain_orientation_file; + grain_orientation_file.push_back(grain_orientation_file_tmp); + Orientation orientation(id, grain_orientation_file, false, inputs.rng_seed, inputs.irf.num_phases, + irf.solidificationTransformation()); + + // Cells at Z = 0 are Solid and do not change type. Cells at Z = 1 and 2 melt at a time step corresponding to their + // Y location in the overall domain (depends on y_offset of the rank) and cool below the liquidus 2 time steps after + // melting. + temperature.num_liquidus_times_this_layer = 2 * grid.nx * grid.ny_local * (grid.nz_layer - 1); + std::vector liquidus_time_step_read(temperature.num_liquidus_times_this_layer), + cell_read(temperature.num_liquidus_times_this_layer); + std::vector cooling_rate_read(temperature.num_liquidus_times_this_layer); + int liq_event_count = 0; + for (int index = 0; index < grid.domain_size; index++) { + grain_id(index) = 1; + celldata.cell_type(index) = Solid; + // Cell coordinates on this rank in X, Y, and Z (GlobalZ = relative to domain bottom) + int coord_z = grid.getCoordZ(index); + int coord_y = grid.getCoordY(index); + // Cells at Z = 0 are Solid, Z = 1 and 2 are Solid but have tm,tl,cr data + if (coord_z > 0) { + // Cells "melt" at a time step corresponding to their Y location in the overall domain (depends on + // y_offset of the rank) + liquidus_time_step_read[liq_event_count] = coord_y + grid.y_offset + 1; + cell_read[liq_event_count] = index; + cooling_rate_read[liq_event_count] = -1.0; + liq_event_count++; + // Cells reach liquidus during cooling 2 time steps after melting + liquidus_time_step_read[liq_event_count] = liquidus_time_step_read[liq_event_count - 1] + 2; + cell_read[liq_event_count] = index; + cooling_rate_read[liq_event_count] = 0.2; + liq_event_count++; + } + } + // Sort event lists and copy to temperature struct + temperature.num_liquidus_times_this_layer = liq_event_count; + temperature.initOrderedTimeTempHistory(liquidus_time_step_read, cell_read, cooling_rate_read, liq_event_count); // Interface struct - int id; - MPI_Comm_rank(MPI_COMM_WORLD, &id); Interface interface(id, grid.domain_size, 0.01); int numcycles = 15; for (int cycle = 1; cycle <= numcycles; cycle++) { - // Update cell types, local undercooling each time step, and fill the steering vector - fillSteeringVector_Remelt(cycle, grid, irf, celldata, temperature, interface); + // Update cell types and fill the steering vector for melting and activation cell type transformations + remeltActivateCells(cycle, grid, irf, celldata, temperature, interface); } // Copy data back to host to check steering vector construction results auto cell_type_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), celldata.cell_type); + Kokkos::deep_copy(interface.num_steer_host, interface.num_steer); auto steering_vector_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), interface.steering_vector); auto num_steer_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), interface.num_steer); auto undercooling_current_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.undercooling_current); - auto liquidus_time_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.liquidus_time); - - // Check the modified cell_type and undercooling_current values on the host: - // Check that the cells corresponding to outside of the "active" portion of the domain have unchanged values - // Check that the cells corresponding to the "active" portion of the domain have potentially changed values - // Z = 3: Rank 0, with all cells having GrainID = 0, should have Liquid cells everywhere, with undercooling - // depending on the Y position Z = 3, Rank > 0, with GrainID > 0, should have TempSolid cells (if not enough time - // steps to melt), Liquid cells (if enough time steps to melt but not reach the liquidus again), or FutureActive - // cells (if below liquidus time step). If the cell is FutureActive type, it should also have a local undercooling - // based on the Rank ID and the time step that it reached the liquidus relative to numcycles Z = 4, all ranks: - // should have TempSolid cells (if not enough time steps to melt) or Liquid (if enough time steps to melt). Local - // undercooling should be based on the Rank ID and the time step that it reached the liquidus relative to numcycles - int future_active_cells = 0; + auto last_time_below_liquidus_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.last_time_below_liquidus); + auto current_cooling_rate_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.current_cooling_rate); + + // Cells at Z = 0 should still be Solid. Z = 1 and 2 cells should be Liquid if numcycles >= the time step the cell + // went above the liquidus, i.e. global_y_cell_position+1, unless the cell is at Z = 2 and if numcycles >= the time + // step the cell cooled below the liquidus, i.e. global_y_cell_position+3, in which case it should be FutureActive + // FutureActive cells should have values for last_time_below_liquidus and current_cooling_rate set accordingly, + // while FutureLiquid cells should have last_time_below_liquidus values of std::numeric_limits::max() for (int index = 0; index < grid.domain_size; index++) { int coord_z = grid.getCoordZ(index); - if (coord_z <= 2) - EXPECT_FLOAT_EQ(undercooling_current_host(index), 0.0); - else { - if (numcycles < liquidus_time_host(index, 0, 0)) { - EXPECT_EQ(cell_type_host(index), TempSolid); - EXPECT_FLOAT_EQ(undercooling_current_host(index), 0.0); - } - else if ((numcycles >= liquidus_time_host(index, 0, 0)) && (numcycles <= liquidus_time_host(index, 0, 1))) { - EXPECT_EQ(cell_type_host(index), Liquid); - EXPECT_FLOAT_EQ(undercooling_current_host(index), 0.0); - } + int coord_y = grid.getCoordY(index); + int expected_liq_time_step_melt = coord_y + grid.y_offset + 1; + int expected_liq_time_step_sol = expected_liq_time_step_melt + 2; + if (coord_z == 0) + EXPECT_EQ(cell_type_host(index), Solid); + else if (coord_z == 1) { + if (numcycles < expected_liq_time_step_melt) + EXPECT_EQ(cell_type_host(index), Solid); else { - EXPECT_FLOAT_EQ(undercooling_current_host(index), (numcycles - liquidus_time_host(index, 0, 1)) * 0.2); - if (coord_z == 4) - EXPECT_EQ(cell_type_host(index), Liquid); - else { + if (numcycles >= expected_liq_time_step_sol) EXPECT_EQ(cell_type_host(index), FutureActive); - future_active_cells++; - } + else + EXPECT_EQ(cell_type_host(index), Liquid); } } + else if (coord_z == 2) { + if (numcycles >= expected_liq_time_step_melt) { + EXPECT_EQ(cell_type_host(index), Liquid); + } + else + EXPECT_EQ(cell_type_host(index), Solid); + } + if (cell_type_host(index) == Liquid) { + if (numcycles >= expected_liq_time_step_sol) + EXPECT_EQ(last_time_below_liquidus_host(index), expected_liq_time_step_sol); + else + EXPECT_EQ(last_time_below_liquidus_host(index), std::numeric_limits::max()); + } + else if (cell_type_host(index) == FutureActive) { + int coord_y = grid.getCoordY(index); + EXPECT_EQ(last_time_below_liquidus_host(index), coord_y + grid.y_offset + 3); + EXPECT_FLOAT_EQ(current_cooling_rate_host(index), 0.2); + } } - // Check the steering vector values on the host - EXPECT_EQ(future_active_cells, num_steer_host(0)); - for (int i = 0; i < future_active_cells; i++) { - // This cell should correspond to a cell at GlobalZ = 3 (RankZ = 1), and some X and Y - int lower_bound_cell_location = grid.nx * grid.ny_local - 1; - int upper_bound_cell_location = 2 * grid.nx * grid.ny_local; - EXPECT_GT(steering_vector_host(i), lower_bound_cell_location); - EXPECT_LT(steering_vector_host(i), upper_bound_cell_location); + // Check that each cell appears on the steering vector the correct number of times: Solid cells = 0 times, + // FutureActive cells = 1 time + std::vector steering_vector_appearences_count(grid.domain_size); + for (int i = 0; i < interface.num_steer_host(0); i++) { + const int index = steering_vector_host(i); + steering_vector_appearences_count[index]++; + } + for (int index = 0; index < grid.domain_size; index++) { + if (cell_type_host(index) == Solid) { + EXPECT_EQ(steering_vector_appearences_count[index], 0); + } + else if (cell_type_host(index) == FutureActive) { + EXPECT_EQ(steering_vector_appearences_count[index], 1); + } } } @@ -649,9 +680,11 @@ void testCalcCritDiagonalLength() { grain_unit_vector = Kokkos::create_mirror_view_and_copy(memory_space(), grain_unit_vector_host); // Initialize interface struct + Grid grid; + grid.domain_size = domain_size; int id; MPI_Comm_rank(MPI_COMM_WORLD, &id); - Interface interface(id, domain_size, 0.01); + Interface interface(id, grid.domain_size, 0.01); // Load octahedron centers into test view view_type octahedron_center_test(Kokkos::ViewAllocateWithoutInitializing("DOCenter"), 3 * domain_size); @@ -742,7 +775,7 @@ TEST(TEST_CATEGORY, communication) { testResizeBuffers(); } TEST(TEST_CATEGORY, cell_update_tests) { - testFillSteeringVector_Remelt(); + testRemeltActivateCells(); testCalcCritDiagonalLength(); testCreateNewOctahedron(); } diff --git a/unit_test/tstNucleation.hpp b/unit_test/tstNucleation.hpp index 83d43f37..4261491e 100644 --- a/unit_test/tstNucleation.hpp +++ b/unit_test/tstNucleation.hpp @@ -38,8 +38,11 @@ void testNucleiInit() { MPI_Comm_rank(MPI_COMM_WORLD, &id); // default inputs struct Inputs inputs; + inputs.domain.deltat = 1.0; // manually set grid for test of 2 layer problem int number_of_layers_temp = 2; + // Initializing layer 1 + int layernumber = 1; Grid grid(number_of_layers_temp); // Create test data // Only Z = 1 through 4 is part of this layer @@ -59,6 +62,7 @@ void testNucleiInit() { grid.y_max = grid.ny * grid.deltax; grid.z_min_layer[1] = grid.z_layer_bottom * grid.deltax; grid.z_max_layer[1] = (grid.z_layer_bottom + grid.nz_layer - 1) * grid.deltax; + grid.layer_height = 1; grid.domain_size = grid.nx * grid.ny_local * grid.nz_layer; grid.domain_size_all_layers = grid.nx * grid.ny_local * grid.nz; grid.bottom_of_current_layer = grid.getBottomOfCurrentLayer(); @@ -75,7 +79,6 @@ void testNucleiInit() { grid.at_north_boundary = false; // Simple IRF - single phase - inputs.domain.deltat = 1 * pow(10, -6); inputs.irf.num_phases = 1; inputs.irf.A[0] = 0.0; inputs.irf.B[0] = 0.0; @@ -97,62 +100,56 @@ void testNucleiInit() { inputs.nucleation.n_max = 0.125; // Allocate temperature data structures - Temperature temperature(grid, inputs.temperature); - // Resize liquidus_time and cooling_rate with the known max number of solidification events - Kokkos::resize(temperature.liquidus_time, grid.domain_size, max_solidification_events_count, 2); - Kokkos::resize(temperature.cooling_rate, grid.domain_size, max_solidification_events_count); - // Initialize max_solidification_events to 3 for each layer. liquidus_time and - // number_of_solidification_events are initialized for each cell on the host and copied to the device - Kokkos::View max_solidification_events_host( - Kokkos::ViewAllocateWithoutInitializing("max_solidification_events_host"), grid.number_of_layers); - max_solidification_events_host(0) = max_solidification_events_count; - max_solidification_events_host(1) = max_solidification_events_count; - // Cells solidify 1, 2, or 3 times, depending on their X coordinate - Kokkos::parallel_for( - "NumSolidificationEventsInit", grid.nz_layer, KOKKOS_LAMBDA(const int &coord_z) { - for (int coord_x = 0; coord_x < grid.nx; coord_x++) { - for (int coord_y = 0; coord_y < grid.ny_local; coord_y++) { - int index = grid.get1DIndex(coord_x, coord_y, coord_z); - if (coord_x < grid.nx / 2 - 1) - temperature.number_of_solidification_events(index) = 3; - else if (coord_x < grid.nx / 2) - temperature.number_of_solidification_events(index) = 2; - else - temperature.number_of_solidification_events(index) = 1; - } + Temperature temperature(grid, inputs.temperature, inputs.print); + // Init first_value, last_value, and raw_temperature_data + // Assume 10 events already happened in layer 0 + temperature.first_value(0) = 0; + temperature.last_value(0) = 10; + int num_temperature_data_pts = 10; + for (int n = 0; n < 3; n++) { + for (int index = 0; index < grid.domain_size; index++) { + // Cell coordinates on this rank in X, Y, and Z (GlobalZ = relative to domain bottom) + int coord_z = grid.getCoordZ(index); + int coord_z_all_layers = coord_z + grid.z_layer_bottom; + int coord_y = grid.getCoordY(index); + int coord_y_global = coord_y + grid.y_offset; + int coord_x = grid.getCoordX(index); + int num_solidification_events; + if (coord_x < grid.nx / 2 - 1) + num_solidification_events = 3; + else if (coord_x < grid.nx / 2) + num_solidification_events = 2; + else + num_solidification_events = 1; + if (n < num_solidification_events) { + temperature.raw_temperature_data(num_temperature_data_pts, 0) = coord_x * grid.deltax; + temperature.raw_temperature_data(num_temperature_data_pts, 1) = coord_y_global * grid.deltax; + temperature.raw_temperature_data(num_temperature_data_pts, 2) = coord_z * grid.deltax; + // melting time step depends on solidification event number + int melt_time = coord_z_all_layers + coord_y + grid.y_offset + (grid.domain_size * n); + temperature.raw_temperature_data(num_temperature_data_pts, 3) = melt_time * inputs.domain.deltat; + // liquidus time stemp depends on solidification event number + temperature.raw_temperature_data(num_temperature_data_pts, 4) = (melt_time + 1) * inputs.domain.deltat; + // ensures that a cell's nucleation time will be 1 time step after its CritTimeStep value + temperature.raw_temperature_data(num_temperature_data_pts, 5) = 1.2; + num_temperature_data_pts++; } - }); - Kokkos::fence(); - Kokkos::parallel_for( - "layer_time_temp_historyInit", max_solidification_events_count, KOKKOS_LAMBDA(const int &n) { - for (int coord_z = 0; coord_z < grid.nz_layer; coord_z++) { - for (int coord_x = 0; coord_x < grid.nx; coord_x++) { - for (int coord_y = 0; coord_y < grid.ny_local; coord_y++) { - int index = grid.get1DIndex(coord_x, coord_y, coord_z); - int coord_z_all_layers = coord_z + grid.z_layer_bottom; - if (n < temperature.number_of_solidification_events(index)) { - // melting time step depends on solidification event number - temperature.liquidus_time(index, n, 0) = - coord_z_all_layers + coord_y + grid.y_offset + (grid.domain_size * n); - // liquidus time stemp depends on solidification event number - temperature.liquidus_time(index, n, 1) = - coord_z_all_layers + coord_y + grid.y_offset + 1 + (grid.domain_size * n); - // ensures that a cell's nucleation time will be 1 time step after its CritTimeStep value - temperature.cooling_rate(index, n) = 1.2; - } - } - } - } - }); - Kokkos::fence(); - temperature.max_solidification_events = - Kokkos::create_mirror_view_and_copy(memory_space(), max_solidification_events_host); + } + } + // Resize with size now known + Kokkos::resize(temperature.raw_temperature_data, num_temperature_data_pts, 6); + temperature.first_value(layernumber) = 10; + temperature.last_value(layernumber) = num_temperature_data_pts; + + // Init temperature data structures (freezing_range value not used in any calculations + double freezing_range = 20.0; + temperature.initialize(layernumber, id, grid, freezing_range, inputs.domain.deltat); // Nucleation data structure, containing views of nuclei locations, time steps, and ids, and nucleation event // counters - initialized with an estimate on the number of nuclei in the layer // (estimated_nuclei_this_rank_this_layer) const double domain_volume = (grid.x_max - grid.x_min) * (grid.y_max - grid.y_min) * - (grid.z_max_layer(1) - grid.z_min_layer(1)) * pow(grid.deltax, 3); + (grid.z_max_layer(layernumber) - grid.z_min_layer(layernumber)) * pow(grid.deltax, 3); int estimated_nuclei_this_rank_this_layer = inputs.nucleation.n_max * pow(grid.deltax, 3) * grid.domain_size; Nucleation nucleation( estimated_nuclei_this_rank_this_layer, inputs.nucleation, @@ -165,7 +162,7 @@ void testNucleiInit() { // Fill in nucleation data structures, and assign nucleation undercooling values to potential nucleation events // Potential nucleation grains are only associated with liquid cells in layer 1 - they will be initialized for each // successive layer when layer 1 in complete - nucleation.placeNuclei(temperature, irf, inputs.rng_seed, 1, grid, id); + nucleation.placeNuclei("FromFile", temperature, irf, inputs.rng_seed, 1, grid, id, inputs.domain.deltat); // Copy results back to host to check auto nuclei_location_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), nucleation.nuclei_locations); @@ -190,10 +187,14 @@ void testNucleiInit() { int coord_y = grid.getCoordY(index); int coord_z_all_layers = coord_z + grid.z_layer_bottom; // Expected nucleation time with remelting can be one of 3 possibilities, depending on the associated - // solidification event - int expected_nucleation_time_no_rm = coord_z_all_layers + coord_y + grid.y_offset + 2; + // solidification event. Should be 1 time step after the liquidus time which in turn is 1 time step after the + // melting time + int melt_time = coord_z_all_layers + coord_y + grid.y_offset; + int expected_nucleation_time_no_rm = melt_time + 2; int associated_s_event = nucleation.nucleation_times_host(n) / grid.domain_size; - int expected_nucleation_time_rm = expected_nucleation_time_no_rm + associated_s_event * grid.domain_size; + // Account for the fact that the melting/liquidus data from the input are increased by 1 time step to avoid + // nucleation prior to the first time step + int expected_nucleation_time_rm = expected_nucleation_time_no_rm + associated_s_event * grid.domain_size + 1; EXPECT_EQ(nucleation.nucleation_times_host(n), expected_nucleation_time_rm); // Are the nucleation events in order of the time steps at which they may occur? diff --git a/unit_test/tstTemperature.hpp b/unit_test/tstTemperature.hpp index 834fe749..0261b7fd 100644 --- a/unit_test/tstTemperature.hpp +++ b/unit_test/tstTemperature.hpp @@ -127,7 +127,7 @@ void testReadTemperatureData(int number_of_layers, bool layerwise_temp_read, boo inputs.temperature.layerwise_temp_read = layerwise_temp_read; // Ensure that constructor correctly initialized the local values of inputs - Temperature temperature(grid, inputs.temperature); + Temperature temperature(grid, inputs.temperature, inputs.print); if (layerwise_temp_read) EXPECT_TRUE(temperature._inputs.layerwise_temp_read); else @@ -222,7 +222,7 @@ void testInit_UnidirectionalGradient(const std::string simulation_type, const do double R_norm = inputs.temperature.R * deltat; // Temperature struct - Temperature temperature(grid, inputs.temperature); + Temperature temperature(grid, inputs.temperature, inputs.print); // Test constructor initialization of _inputs // These should've been initialized with default values EXPECT_FALSE(temperature._inputs.layerwise_temp_read); @@ -234,16 +234,14 @@ void testInit_UnidirectionalGradient(const std::string simulation_type, const do temperature.initialize(id, simulation_type, grid, deltat); // Copy temperature views back to host - auto number_of_solidification_events_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.number_of_solidification_events); - auto solidification_event_counter_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.solidification_event_counter); - auto liquidus_time_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.liquidus_time); - auto max_solidification_events_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.max_solidification_events); - auto undercooling_current_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.undercooling_current); - auto cooling_rate_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.cooling_rate); + auto liquidus_cell_list_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.liquidus_cell_list); + auto cooling_rate_list_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.cooling_rate_list); + auto last_time_below_liquidus_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.last_time_below_liquidus); + auto current_cooling_rate_host = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), temperature.current_cooling_rate); // Check results int location_init_undercooling, location_liquidus_isotherm; @@ -260,42 +258,38 @@ void testInit_UnidirectionalGradient(const std::string simulation_type, const do location_liquidus_isotherm = location_init_undercooling + Kokkos::round(inputs.temperature.init_undercooling / (G * grid.deltax)); - EXPECT_EQ(max_solidification_events_host(0), 1); + // Each cell melts/solidifies once + EXPECT_EQ(temperature.num_liquidus_times_this_layer, 2 * grid.domain_size); + // The counter is at num_liquidus_times_this_layer since melting/resolidification aren't considered for this problem + // type + EXPECT_EQ(temperature.num_liquidus_times_this_layer, temperature.liquidus_time_counter); + std::vector found_liq_events_per_cell(grid.domain_size, 0); + for (int n = 0; n < temperature.num_liquidus_times_this_layer; n++) { + // Each cell should appear on this list twice, once for melting (time step 0, cooling rate -1) and once for + // cooling below the liquidus (cooling rate = 1, liquidus time either 0 if G = 0 or a Z dependent value + const int index = liquidus_cell_list_host(n); + found_liq_events_per_cell[index]++; + if (cooling_rate_list_host[n] == -1.0) { + EXPECT_EQ(temperature.liquidus_time_list_host[n], 0); + } + else { + EXPECT_FLOAT_EQ(cooling_rate_list_host[n], 1.0); + const int coord_z = grid.getCoordZ(index); + const int expected_liq_time = Kokkos::round((coord_z - location_liquidus_isotherm) * G_norm / R_norm); + EXPECT_EQ(temperature.liquidus_time_list_host(n), expected_liq_time); + } + } + // Check correct initialization of current cooling rate and last time below liquidus, which are the views used for + // this problem type to update the temperature field for (int coord_z = 0; coord_z < grid.nz; coord_z++) { for (int coord_x = 0; coord_x < grid.nx; coord_x++) { for (int coord_y = 0; coord_y < grid.ny_local; coord_y++) { int index = grid.get1DIndex(coord_x, coord_y, coord_z); - // Each cell solidifies once, and counter should start at 0, associated with the zeroth layer - // melt time step should be -1 for all cells - // Cells cool at 1 K per time step - EXPECT_FLOAT_EQ(liquidus_time_host(index, 0, 0), -1.0); - EXPECT_FLOAT_EQ(cooling_rate_host(index, 0), R_norm); - EXPECT_EQ(number_of_solidification_events_host(index), 1); - EXPECT_EQ(solidification_event_counter_host(index), 0); - // undercooling_current should be zero for cells if a positive G is given, or init_undercooling if being - // initialized with a uniform undercooling field (all cells initially below liquidus) - if (G == 0) { - EXPECT_FLOAT_EQ(undercooling_current_host(index), inputs.temperature.init_undercooling); - EXPECT_FLOAT_EQ(liquidus_time_host(index, 0, 1), -1); - } - else { - int dist_from_liquidus = coord_z - location_liquidus_isotherm; - if (dist_from_liquidus < 0) { - // Undercooled cell (liquidus time step already passed, set to -1) - // The undercooling at Z = locationOfLiquidus should have been set to init_undercooling - EXPECT_FLOAT_EQ(liquidus_time_host(index, 0, 1), -1); - int dist_from_init_undercooling = coord_z - location_init_undercooling; - EXPECT_FLOAT_EQ(undercooling_current_host(index), - inputs.temperature.init_undercooling - - dist_from_init_undercooling * G * grid.deltax); - } - else { - // Cell has not yet reached the nonzero liquidus time yet (or reaches it at time = 0), either - // does not have an assigned undercooling or the assigned undercooling is zero - EXPECT_FLOAT_EQ(liquidus_time_host(index, 0, 1), dist_from_liquidus * G_norm / R_norm); - EXPECT_FLOAT_EQ(undercooling_current_host(index), 0.0); - } - } + const int expected_liq_time = Kokkos::round((coord_z - location_liquidus_isotherm) * G_norm / R_norm); + EXPECT_EQ(last_time_below_liquidus_host(index), expected_liq_time); + EXPECT_FLOAT_EQ(current_cooling_rate_host(index), R_norm); + EXPECT_EQ(temperature.number_of_solidification_events(index), 1); + EXPECT_EQ(found_liq_events_per_cell[index], 2); } } } @@ -428,7 +422,7 @@ void testInitTemperatureFromFinch(const bool mirror_x, const int number_of_copie } // Construct temperature views - Temperature temperature(id, np, grid, inputs.temperature, input_temperature_data); + Temperature temperature(id, np, grid, inputs.temperature, inputs.print, input_temperature_data); // Check that the right ranks have the right temperature data checkTemperatureResults(inputs, grid, temperature, number_of_copies, id, np, mirror_x); @@ -477,7 +471,7 @@ void testInitTemperatureFromFile(const bool mirror_x, const int number_of_copies inputs.temperature.temp_files_in_series = 1; inputs.temperature.layerwise_temp_read = false; - Temperature temperature(grid, inputs.temperature); + Temperature temperature(grid, inputs.temperature, inputs.print); temperature.readTemperatureData(id, grid, 0); // Check that the right ranks have the right temperature data diff --git a/unit_test/tstUpdate.hpp b/unit_test/tstUpdate.hpp index 181f8dad..de13b4ae 100644 --- a/unit_test/tstUpdate.hpp +++ b/unit_test/tstUpdate.hpp @@ -43,7 +43,7 @@ void testSmallDirS() { Grid grid(inputs.simulation_type, id, np, inputs.domain.number_of_layers, inputs.domain, inputs.substrate, inputs.temperature); // Temperature fields characterized by data in this structure - Temperature temperature(grid, inputs.temperature, inputs.print.store_solidification_start); + Temperature temperature(grid, inputs.temperature, inputs.print); // Run SmallDirS problem and check volume fraction of nucleated grains with 1% tolerance of expected value (to // account for the non-deterministic nature of the cell capture) @@ -75,7 +75,7 @@ void testSmallEquiaxedGrain() { Grid grid(inputs.simulation_type, id, np, inputs.domain.number_of_layers, inputs.domain, inputs.substrate, inputs.temperature); // Temperature fields characterized by data in this structure - Temperature temperature(grid, inputs.temperature, inputs.print.store_solidification_start); + Temperature temperature(grid, inputs.temperature, inputs.print); // Run Small equiaxed grain problem and check time step at which the grain reaches the domain edge runExaCA(id, np, inputs, timers, grid, temperature); @@ -86,7 +86,8 @@ void testSmallEquiaxedGrain() { std::ifstream log_data_stream(log_file); nlohmann::json log_data = nlohmann::json::parse(log_data_stream); int time_step_of_output = log_data["TimeStepOfOutput"]; - // FIXME: Output time step is usually 4820, but may be 4821 - need to investigate this possible race condition + // FIXME: Output time step is usually 4820, but may be 1 time step larger - need to investigate this possible race + // condition. Finishes 1 time step earlier now after update to active cell init for no remelting case EXPECT_NEAR(time_step_of_output, 4820, 1); } //---------------------------------------------------------------------------//