diff --git a/libnufi_lib.a b/libnufi_lib.a index 031c6a7..e0677ec 100644 Binary files a/libnufi_lib.a and b/libnufi_lib.a differ diff --git a/nufi/grids.h b/nufi/grids.h index 66a5e0b..0c97fc5 100644 --- a/nufi/grids.h +++ b/nufi/grids.h @@ -2,11 +2,13 @@ #define GRIDS_H #include "nufi/cells.h" +#include #include #include #include #include #include +#include #include #include @@ -25,7 +27,7 @@ template struct GridStructure { CellLocator locator; - unsigned int grid_version; + unsigned int grid_version = 0; // === // === // // Evaluator // @@ -35,49 +37,68 @@ template struct GridStructure { const std::vector> &points) const { std::vector values(points.size()); -#pragma omp parallel - { - std::vector local_solution_buffer(fe->n_dofs_per_cell()); - FEPointEvaluation evaluator(*mapping, *fe, update_gradients); -#pragma omp for - for (unsigned int p = 0; p < points.size(); ++p) { - const auto cell_location = locator.locate(points[p]); +#pragma omp parallel for + for (unsigned int p = 0; p < points.size(); ++p) { - cell_location.info->cell->get_dof_values(solution, - local_solution_buffer.begin(), - local_solution_buffer.end()); + const Tensor<1, dim> grad_phi = VectorTools::point_gradient( + *mapping, *dof_handler, solution, points[p]); - evaluator.reinit( - cell_location.info->cell, - ArrayView>(&cell_location.reference_point, 1)); - - evaluator.evaluate(local_solution_buffer, EvaluationFlags::gradients); - - values[p] = evaluator.get_gradient(0)[0]; - } + values[p] = grad_phi[0]; } + return values; } + // #pragma omp parallel + // { + // std::vector local_solution_buffer(fe->n_dofs_per_cell()); + // FEPointEvaluation evaluator(*mapping, *fe, + // update_gradients); + // #pragma omp for + // for (unsigned int p = 0; p < points.size(); ++p) { + // + // const auto cell_location = locator.locate(points[p]); + // + // cell_location.info->cell->get_dof_values(solution, + // local_solution_buffer.begin(), + // local_solution_buffer.end()); + // + // evaluator.reinit( + // cell_location.info->cell, + // ArrayView>(&cell_location.reference_point, + // 1)); + // + // evaluator.evaluate(local_solution_buffer, + // EvaluationFlags::gradients); + // + // values[p] = evaluator.get_gradient(0)[0]; + // } + // } + // AssertThrow(dof_handler->n_dofs() == solution.size(), + // ExcMessage("Solution vector size does not match + // DoFHandler.")); + // return values; + // } }; template GridStructure make_grid_snapshot(const PoissonProblem &poisson) { GridStructure grid; - grid.grid_version = 0; grid.triangulation = std::make_unique>(); + grid.triangulation->copy_triangulation(poisson.get_triangulation()); + grid.fe = std::make_unique>(poisson.get_fe()); + + grid.dof_handler = std::make_unique>(*grid.triangulation); + + grid.dof_handler->distribute_dofs(*grid.fe); + grid.mapping = std::make_unique>(poisson.get_mapping()); - grid.dof_handler = std::make_unique>(*grid.triangulation); - grid.dof_handler->distribute_dofs(poisson.get_dof_handler().get_fe()); - grid.locator.rebuild(*grid.dof_handler, *grid.triangulation); - grid.fe = std::make_unique>(poisson.get_fe()); - return grid; } diff --git a/nufi/parameters.h b/nufi/parameters.h index 4bcf0be..82e0d5d 100644 --- a/nufi/parameters.h +++ b/nufi/parameters.h @@ -35,10 +35,10 @@ constexpr double F0_FACTOR = 0.39894228040143267793994; // 1/sqrt(2pi) // NUFI options constexpr double DT = 1. / 10.; constexpr unsigned int TMAX = 100; -constexpr unsigned int REFINE_FREQUENCY = 10; +constexpr unsigned int REFINE_FREQUENCY = 3; // Plotting options -constexpr int PLOT_FREQUENCY = 10; +constexpr int PLOT_FREQUENCY = 1; constexpr size_t PLOT_NX = CALC_NX; constexpr double PLOT_DX = LX / PLOT_NX; const std::string PLOT_DIR = "results/"; diff --git a/src/nufi_solver.cc b/src/nufi_solver.cc index be03cfb..60cf187 100644 --- a/src/nufi_solver.cc +++ b/src/nufi_solver.cc @@ -191,8 +191,8 @@ void NuFISolver::run() { std::vector> phi_history; update_grid_versions(grid_versions, poisson); - update_solution_history(phi_history, poisson, - grid_versions.back().grid_version); + // update_solution_history(phi_history, poisson, + // grid_versions.back().grid_version); std::vector x_eval(Parameters::CALC_NX); @@ -239,7 +239,21 @@ void NuFISolver::run() { // } // // std::cout << "h ratio = " << max_h / min_h << std::endl; - + // std::vector x = make_x_eval(Parameters::CALC_NX); + // auto rho = eval_rho(it, x, grid_versions, phi_history, Parameters::NV); + // + // double mean = 0; + // + // for (auto r : rho) + // mean += r; + // + // mean /= rho.size(); + // + // std::cout << "rho mean = " << mean << "\n"; + // std::cout << "rho min = " << *std::min_element(rho.begin(), rho.end()) + // << "\n"; + // std::cout << "rho max = " << *std::max_element(rho.begin(), rho.end()) + // << "\n"; // END: diagnostics double compute_start = timer.elapsed(); @@ -277,24 +291,27 @@ void NuFISolver::run() { double step_time = timer_elapsed - time_elapsed_before; std::cout << "step made in " << step_time << " seconds\n\n"; + + //====//====// + // Plotting // + //====//====// if (it % Parameters::PLOT_FREQUENCY == 0) { double plot_start = timer.elapsed(); + std::cout << "Saving results... "; save_f(*this, it, grid_versions, phi_history, Parameters::PLOT_NX, Parameters::NV, "results/ftilda_" + std::to_string(it) + ".dat"); save_rho(*this, it, grid_versions, phi_history, Parameters::PLOT_NX, "results/rho_" + std::to_string(it) + ".dat"); - // save_Efield(it, coeffs.get(), 128, "results/field_" + - // std::to_string(it) + ".dat"); std::vector x_eval_Ex = make_x_eval(Parameters::PLOT_NX); - std::vector tmp_rho(x_eval_Ex.size()); - tmp_rho = eval(x_eval_Ex, grid_versions[phi_history[it].grid_version], - phi_history[it].solution); + std::vector tmp_Ex(x_eval_Ex.size()); + tmp_Ex = eval(x_eval_Ex, grid_versions[phi_history[it].grid_version], + phi_history[it].solution); std::vector E_x(Parameters::PLOT_NX); for (size_t i = 0; i < Parameters::PLOT_NX; ++i) - E_x[i] = -tmp_rho[i]; + E_x[i] = -tmp_Ex[i]; save_space_vector(E_x, "field", it); double int_val = 0.5 * integral_space_vector_squared(