From 58ca99f9c57cd85fcb8a97bb094b481dd35522db Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 15 Mar 2025 14:42:26 -0400 Subject: [PATCH 1/3] deposit test pgens --- setups/tests/deposit/deposit-gr.toml | 53 +++++++++++ setups/tests/deposit/deposit-mink.toml | 69 +++++++++++++++ setups/tests/deposit/deposit-sr.toml | 53 +++++++++++ setups/tests/deposit/pgen.hpp | 116 +++++++++---------------- 4 files changed, 218 insertions(+), 73 deletions(-) create mode 100644 setups/tests/deposit/deposit-gr.toml create mode 100644 setups/tests/deposit/deposit-mink.toml create mode 100644 setups/tests/deposit/deposit-sr.toml diff --git a/setups/tests/deposit/deposit-gr.toml b/setups/tests/deposit/deposit-gr.toml new file mode 100644 index 000000000..077e3e59a --- /dev/null +++ b/setups/tests/deposit/deposit-gr.toml @@ -0,0 +1,53 @@ +[simulation] + name = "deposit-test" + engine = "srpic" + runtime = 10.0 + +[grid] + resolution = [256, 256] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 0.1 + skindepth0 = 0.1 + +[algorithms] + current_filters = 4 + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 10.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e2 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e2 + +[setup] + +[output] + format = "hdf5" + interval_time = 0.01 + + [output.quantities] + quantities = ["N_1", "N_2", "E", "B", "J"] + +[diagnostics] + colored_stdout = true + blocking_timers = true diff --git a/setups/tests/deposit/deposit-mink.toml b/setups/tests/deposit/deposit-mink.toml new file mode 100644 index 000000000..df70f7552 --- /dev/null +++ b/setups/tests/deposit/deposit-mink.toml @@ -0,0 +1,69 @@ +[simulation] + name = "deposit-test-mink" + engine = "srpic" + runtime = 10.0 + +[grid] + resolution = [512, 512] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 0.1 + skindepth0 = 0.1 + +[algorithms] + current_filters = 4 + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 10.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e2 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e2 + +[setup] + x1s = [0.25] + y1s = [0.85] + z1s = [0.33] + ux1s = [0.6] + uy1s = [-0.3] + uz1s = [-0.2] + + x2s = [0.25] + y2s = [0.85] + z2s = [0.33] + ux2s = [-0.2] + uy2s = [-0.2] + uz2s = [0.1] + +[output] + format = "hdf5" + interval_time = 0.01 + + [output.fields] + quantities = ["N_1", "N_2", "E", "B", "J"] + +[checkpoint] + keep = 0 + +[diagnostics] + colored_stdout = true + blocking_timers = true diff --git a/setups/tests/deposit/deposit-sr.toml b/setups/tests/deposit/deposit-sr.toml new file mode 100644 index 000000000..077e3e59a --- /dev/null +++ b/setups/tests/deposit/deposit-sr.toml @@ -0,0 +1,53 @@ +[simulation] + name = "deposit-test" + engine = "srpic" + runtime = 10.0 + +[grid] + resolution = [256, 256] + extent = [[0.0, 1.0], [0.0, 1.0]] + + [grid.metric] + metric = "minkowski" + + [grid.boundaries] + fields = [["PERIODIC"], ["PERIODIC"]] + particles = [["PERIODIC"], ["PERIODIC"]] + +[scales] + larmor0 = 0.1 + skindepth0 = 0.1 + +[algorithms] + current_filters = 4 + + [algorithms.timestep] + CFL = 0.5 + +[particles] + ppc0 = 10.0 + + [[particles.species]] + label = "e-" + mass = 1.0 + charge = -1.0 + maxnpart = 1e2 + + [[particles.species]] + label = "e+" + mass = 1.0 + charge = 1.0 + maxnpart = 1e2 + +[setup] + +[output] + format = "hdf5" + interval_time = 0.01 + + [output.quantities] + quantities = ["N_1", "N_2", "E", "B", "J"] + +[diagnostics] + colored_stdout = true + blocking_timers = true diff --git a/setups/tests/deposit/pgen.hpp b/setups/tests/deposit/pgen.hpp index fd9a41c2e..1be4deb69 100644 --- a/setups/tests/deposit/pgen.hpp +++ b/setups/tests/deposit/pgen.hpp @@ -4,14 +4,8 @@ #include "enums.h" #include "global.h" -#include "arch/kokkos_aliases.h" #include "arch/traits.h" -#include "utils/comparators.h" -#include "utils/formatting.h" -#include "utils/log.h" -#include "utils/numeric.h" -#include "archetypes/energy_dist.h" #include "archetypes/particle_injector.h" #include "archetypes/problem_generator.h" #include "framework/domain/domain.h" @@ -26,10 +20,20 @@ namespace user { struct PGen : public arch::ProblemGenerator { // compatibility traits for the problem generator - static constexpr auto engines = traits::compatible_with::value; - static constexpr auto metrics = traits::compatible_with::value; - static constexpr auto dimensions = - traits::compatible_with::value; + static constexpr auto engines { + traits::compatible_with::value + }; + static constexpr auto metrics { + traits::compatible_with::value + }; + static constexpr auto dimensions { + traits::compatible_with::value + }; // for easy access to variables in the child class using arch::ProblemGenerator::D; @@ -47,85 +51,51 @@ namespace user { const auto x1s = params.template get>("setup.x1s", empty); const auto y1s = params.template get>("setup.y1s", empty); const auto z1s = params.template get>("setup.z1s", empty); - const auto ux1s = params.template get>("setup.ux1s", + const auto phi1s = params.template get>("setup.phi1s", + empty); + const auto ux1s = params.template get>("setup.ux1s", empty); - const auto uy1s = params.template get>("setup.uy1s", + const auto uy1s = params.template get>("setup.uy1s", empty); - const auto uz1s = params.template get>("setup.uz1s", + const auto uz1s = params.template get>("setup.uz1s", empty); const auto x2s = params.template get>("setup.x2s", empty); const auto y2s = params.template get>("setup.y2s", empty); const auto z2s = params.template get>("setup.z2s", empty); - const auto ux2s = params.template get>("setup.ux2s", + const auto ux2s = params.template get>("setup.ux2s", empty); - const auto uy2s = params.template get>("setup.uy2s", + const auto uy2s = params.template get>("setup.uy2s", empty); - const auto uz2s = params.template get>("setup.uz2s", + const auto uz2s = params.template get>("setup.uz2s", empty); - // std::vector x, y, z, ux_1, uy_1, uz_1, ux_2, uy_2, uz_2; - // x.push_back(0.85); - // x.push_back(0.123); - // if constexpr (D == Dim::_2D || D == Dim::_3D) { - // y.push_back(0.32); - // y.push_back(0.321); - // } - // if constexpr (D == Dim::_3D) { - // z.push_back(0.231); - // z.push_back(0.687); - // } - // ux_1.push_back(1.0); - // uy_1.push_back(-1.0); - // uz_1.push_back(0.0); - // ux_1.push_back(1.0); - // uy_1.push_back(-2.0); - // uz_1.push_back(1.0); - // - // ux_2.push_back(1.0); - // uy_2.push_back(1.0); - // uz_2.push_back(0.0); - // ux_2.push_back(-2.0); - // uy_2.push_back(3.0); - // uz_2.push_back(-1.0); - // - const std::map> data_1 { - { "x1", x1s}, - { "x2", y1s}, - { "x3", z1s}, - {"ux1", ux1s}, - {"ux2", uy1s}, - {"ux3", uz1s} + const auto phi2s = params.template get>("setup.phi2s", + empty); + std::map> data_1 { + { "x1", x1s }, + { "x2", y1s }, + { "ux1", ux1s }, + { "ux2", uy1s }, + { "ux3", uz1s } }; - const std::map> data_2 { - { "x1", x2s}, - { "x2", y2s}, - { "x3", z2s}, - {"ux1", ux2s}, - {"ux2", uy2s}, - {"ux3", uz2s} + std::map> data_2 { + { "x1", x2s }, + { "x2", y2s }, + { "ux1", ux2s }, + { "ux2", uy2s }, + { "ux3", uz2s } }; + if constexpr (M::CoordType == Coord::Cart or D == Dim::_3D) { + data_1["x3"] = z1s; + data_2["x3"] = z2s; + } else if constexpr (D == Dim::_2D) { + data_1["phi"] = phi1s; + data_2["phi"] = phi2s; + } arch::InjectGlobally(global_domain, local_domain, (arch::spidx_t)1, data_1); arch::InjectGlobally(global_domain, local_domain, (arch::spidx_t)2, data_2); } - - // void CustomPostStep(std::size_t, long double time, Domain& domain) { - // if (time >= 0.1) { - // for (auto& species : domain.species) { - // auto ux1 = species.ux1; - // auto ux2 = species.ux2; - // auto ux3 = species.ux3; - // Kokkos::parallel_for( - // "Stop", - // species.rangeActiveParticles(), - // Lambda(index_t p) { - // ux1(p) = ZERO; - // ux2(p) = ZERO; - // ux3(p) = ZERO; - // }); - // } - // } - // } }; } // namespace user From 9fcca0b05f4b6378509cc268131a7bdd20c92bbb Mon Sep 17 00:00:00 2001 From: haykh Date: Sat, 15 Mar 2025 14:43:12 -0400 Subject: [PATCH 2/3] output domain shapes --- src/framework/domain/communications.cpp | 18 +++++----- src/framework/domain/output.cpp | 18 +++++++++- src/output/tests/writer-mpi.cpp | 1 + src/output/tests/writer-nompi.cpp | 1 + src/output/writer.cpp | 48 ++++++++++++++----------- src/output/writer.h | 6 +++- 6 files changed, 60 insertions(+), 32 deletions(-) diff --git a/src/framework/domain/communications.cpp b/src/framework/domain/communications.cpp index 7dc5d285a..0bd06b82e 100644 --- a/src/framework/domain/communications.cpp +++ b/src/framework/domain/communications.cpp @@ -36,10 +36,10 @@ namespace ntt { using comm_params_t = std::pair>; template - auto GetSendRecvRanks(Metadomain* metadomain, - Domain& domain, - dir::direction_t direction) - -> std::pair { + auto GetSendRecvRanks( + Metadomain* metadomain, + Domain& domain, + dir::direction_t direction) -> std::pair { Domain* send_to_nghbr_ptr = nullptr; Domain* recv_from_nghbr_ptr = nullptr; // set pointers to the correct send/recv domains @@ -119,11 +119,11 @@ namespace ntt { } template - auto GetSendRecvParams(Metadomain* metadomain, - Domain& domain, - dir::direction_t direction, - bool synchronize) - -> std::pair { + auto GetSendRecvParams( + Metadomain* metadomain, + Domain& domain, + dir::direction_t direction, + bool synchronize) -> std::pair { const auto [send_indrank, recv_indrank] = GetSendRecvRanks(metadomain, domain, direction); const auto [send_ind, send_rank] = send_indrank; diff --git a/src/framework/domain/output.cpp b/src/framework/domain/output.cpp index 3dc95c18a..c51894e55 100644 --- a/src/framework/domain/output.cpp +++ b/src/framework/domain/output.cpp @@ -70,6 +70,7 @@ namespace ntt { g_writer.defineMeshLayout(glob_shape_with_ghosts, off_ncells_with_ghosts, loc_shape_with_ghosts, + { local_domain->index(), ndomains() }, params.template get>( "output.fields.downsampling"), incl_ghosts, @@ -227,6 +228,17 @@ namespace ntt { const auto dwn = params.template get>( "output.fields.downsampling"); + auto off_ncells_with_ghosts = local_domain->offset_ncells(); + auto loc_shape_with_ghosts = local_domain->mesh.n_active(); + { // compute positions/sizes of meshblocks in cells in all dimensions + const auto off_ndomains = local_domain->offset_ndomains(); + if (incl_ghosts) { + for (auto d { 0 }; d <= M::Dim; ++d) { + off_ncells_with_ghosts[d] += 2 * N_GHOSTS * off_ndomains[d]; + loc_shape_with_ghosts[d] += 2 * N_GHOSTS; + } + } + } for (unsigned short dim = 0; dim < M::Dim; ++dim) { const auto l_size = local_domain->mesh.n_active()[dim]; const auto l_offset = local_domain->offset_ncells()[dim]; @@ -275,7 +287,11 @@ namespace ntt { xe(offset + i_dwn + 1) = x_Ph[dim]; } }); - g_writer.writeMesh(dim, xc, xe); + g_writer.writeMesh( + dim, + xc, + xe, + { off_ncells_with_ghosts[dim], loc_shape_with_ghosts[dim] }); } const auto output_asis = params.template get("output.debug.as_is"); // !TODO: this can probably be optimized to dump things at once diff --git a/src/output/tests/writer-mpi.cpp b/src/output/tests/writer-mpi.cpp index f6d3ee88a..4ec6941ac 100644 --- a/src/output/tests/writer-mpi.cpp +++ b/src/output/tests/writer-mpi.cpp @@ -64,6 +64,7 @@ auto main(int argc, char* argv[]) -> int { writer.defineMeshLayout({ static_cast(mpi_size) * nx1 }, { static_cast(mpi_rank) * nx1 }, { nx1 }, + { mpi_rank, mpi_size }, { dwn1 }, false, Coord::Cart); diff --git a/src/output/tests/writer-nompi.cpp b/src/output/tests/writer-nompi.cpp index 8fb2ac026..7e795c2a9 100644 --- a/src/output/tests/writer-nompi.cpp +++ b/src/output/tests/writer-nompi.cpp @@ -74,6 +74,7 @@ auto main(int argc, char* argv[]) -> int { writer.defineMeshLayout({ nx1, nx2, nx3 }, { 0, 0, 0 }, { nx1, nx2, nx3 }, + { 0, 1 }, { dwn1, dwn2, dwn3 }, false, Coord::Cart); diff --git a/src/output/writer.cpp b/src/output/writer.cpp index 95965c864..e2e094535 100644 --- a/src/output/writer.cpp +++ b/src/output/writer.cpp @@ -63,12 +63,14 @@ namespace out { m_mode = mode; } - void Writer::defineMeshLayout(const std::vector& glob_shape, - const std::vector& loc_corner, - const std::vector& loc_shape, - const std::vector& dwn, - bool incl_ghosts, - Coord coords) { + void Writer::defineMeshLayout( + const std::vector& glob_shape, + const std::vector& loc_corner, + const std::vector& loc_shape, + const std::pair& domain_idx, + const std::vector& dwn, + bool incl_ghosts, + Coord coords) { m_flds_ghosts = incl_ghosts; m_dwn = dwn; @@ -99,24 +101,24 @@ namespace out { for (std::size_t i { 0 }; i < m_flds_g_shape.size(); ++i) { // cell-centers - adios2::Dims g_shape = { m_flds_g_shape_dwn[i] }; - adios2::Dims l_corner = { m_flds_l_corner_dwn[i] }; - adios2::Dims l_shape = { m_flds_l_shape_dwn[i] }; m_io.DefineVariable("X" + std::to_string(i + 1), - g_shape, - l_corner, - l_shape, + { m_flds_g_shape_dwn[i] }, + { m_flds_l_corner_dwn[i] }, + { m_flds_l_shape_dwn[i] }, adios2::ConstantDims); // cell-edges - const auto is_last = (m_flds_l_corner[i] + m_flds_l_shape[i] == + const auto is_last = (m_flds_l_corner[i] + m_flds_l_shape[i] == m_flds_g_shape[i]); - adios2::Dims g_shape1 = { m_flds_g_shape_dwn[i] + 1 }; - adios2::Dims l_shape1 = { m_flds_l_shape_dwn[i] + (is_last ? 1 : 0) }; m_io.DefineVariable("X" + std::to_string(i + 1) + "e", - g_shape1, - l_corner, - l_shape1, + { m_flds_g_shape_dwn[i] + 1 }, + { m_flds_l_corner_dwn[i] }, + { m_flds_l_shape_dwn[i] + (is_last ? 1 : 0) }, adios2::ConstantDims); + m_io.DefineVariable("N" + std::to_string(i + 1) + "l", + { 2 * domain_idx.second }, + { 2 * domain_idx.first }, + { 2 }, + adios2::ConstantDims); } if constexpr (std::is_same::array_layout, @@ -402,9 +404,10 @@ namespace out { m_writer.Put(var, e_bins_h); } - void Writer::writeMesh(unsigned short dim, - const array_t& xc, - const array_t& xe) { + void Writer::writeMesh(unsigned short dim, + const array_t& xc, + const array_t& xe, + const std::vector& loc_off_sz) { auto varc = m_io.InquireVariable("X" + std::to_string(dim + 1)); auto vare = m_io.InquireVariable("X" + std::to_string(dim + 1) + "e"); auto xc_h = Kokkos::create_mirror_view(xc); @@ -413,6 +416,9 @@ namespace out { Kokkos::deep_copy(xe_h, xe); m_writer.Put(varc, xc_h); m_writer.Put(vare, xe_h); + auto vard = m_io.InquireVariable( + "N" + std::to_string(dim + 1) + "l"); + m_writer.Put(vard, loc_off_sz.data()); } void Writer::beginWriting(WriteModeTags write_mode, diff --git a/src/output/writer.h b/src/output/writer.h index a8abf4b12..0b95dfb0e 100644 --- a/src/output/writer.h +++ b/src/output/writer.h @@ -86,6 +86,7 @@ namespace out { void defineMeshLayout(const std::vector&, const std::vector&, const std::vector&, + const std::pair&, const std::vector&, bool, Coord); @@ -94,7 +95,10 @@ namespace out { void defineParticleOutputs(Dimension, const std::vector&); void defineSpectraOutputs(const std::vector&); - void writeMesh(unsigned short, const array_t&, const array_t&); + void writeMesh(unsigned short, + const array_t&, + const array_t&, + const std::vector&); template void writeField(const std::vector&, From 1e18874d3c870198ba1c74ead9ab59c1f0eca51b Mon Sep 17 00:00:00 2001 From: haykh Date: Thu, 27 Mar 2025 17:45:57 -0400 Subject: [PATCH 3/3] i+di coord push for cart --- src/kernels/particle_pusher_sr.hpp | 112 +++++++++++++++++++---------- 1 file changed, 76 insertions(+), 36 deletions(-) diff --git a/src/kernels/particle_pusher_sr.hpp b/src/kernels/particle_pusher_sr.hpp index 2e8a5f652..50fbec5b0 100644 --- a/src/kernels/particle_pusher_sr.hpp +++ b/src/kernels/particle_pusher_sr.hpp @@ -562,45 +562,85 @@ namespace kernel::sr { Inline void posUpd(bool massive, index_t& p, coord_t& xp) const { // get cartesian velocity - const real_t inv_energy { - massive ? ONE / math::sqrt(ONE + SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) - : ONE / math::sqrt(SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) - }; - vec_t vp_Cart { ux1(p) * inv_energy, - ux2(p) * inv_energy, - ux3(p) * inv_energy }; - // get cartesian position - coord_t xp_Cart { ZERO }; - metric.template convert_xyz(xp, xp_Cart); - // update cartesian position - for (auto d = 0u; d < M::PrtlDim; ++d) { - xp_Cart[d] += vp_Cart[d] * dt; - } - // transform back to code - metric.template convert_xyz(xp_Cart, xp); - - // update x1 - if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) { - i1_prev(p) = i1(p); - dx1_prev(p) = dx1(p); - from_Xi_to_i_di(xp[0], i1(p), dx1(p)); - } + if constexpr (M::CoordType == Coord::Cart) { + // i+di push for Cartesian basis + const real_t dt_inv_energy { + massive + ? (dt / math::sqrt(ONE + SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p)))) + : (dt / math::sqrt(SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p)))) + }; + if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) { + i1_prev(p) = i1(p); + dx1_prev(p) = dx1(p); + dx1(p) += metric.template transform<1, Idx::XYZ, Idx::U>(xp, ux1(p)) * + dt_inv_energy; + i1(p) += static_cast(dx1(p) >= ONE) - + static_cast(dx1(p) < ZERO); + dx1(p) -= (dx1(p) >= ONE); + dx1(p) += (dx1(p) < ZERO); + } + if constexpr (D == Dim::_2D || D == Dim::_3D) { + i2_prev(p) = i2(p); + dx2_prev(p) = dx2(p); + dx2(p) += metric.template transform<2, Idx::XYZ, Idx::U>(xp, ux2(p)) * + dt_inv_energy; + i2(p) += static_cast(dx2(p) >= ONE) - + static_cast(dx2(p) < ZERO); + dx2(p) -= (dx2(p) >= ONE); + dx2(p) += (dx2(p) < ZERO); + } + if constexpr (D == Dim::_3D) { + i3_prev(p) = i3(p); + dx3_prev(p) = dx3(p); + dx3(p) += metric.template transform<3, Idx::XYZ, Idx::U>(xp, ux3(p)) * + dt_inv_energy; + i3(p) += static_cast(dx3(p) >= ONE) - + static_cast(dx3(p) < ZERO); + dx3(p) -= (dx3(p) >= ONE); + dx3(p) += (dx3(p) < ZERO); + } + } else { + // full Cartesian coordinate push in non-Cartesian basis + const real_t inv_energy { + massive ? ONE / math::sqrt(ONE + SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) + : ONE / math::sqrt(SQR(ux1(p)) + SQR(ux2(p)) + SQR(ux3(p))) + }; + vec_t vp_Cart { ux1(p) * inv_energy, + ux2(p) * inv_energy, + ux3(p) * inv_energy }; + // get cartesian position + coord_t xp_Cart { ZERO }; + metric.template convert_xyz(xp, xp_Cart); + // update cartesian position + for (auto d = 0u; d < M::PrtlDim; ++d) { + xp_Cart[d] += vp_Cart[d] * dt; + } + // transform back to code + metric.template convert_xyz(xp_Cart, xp); + + // update x1 + if constexpr (D == Dim::_1D || D == Dim::_2D || D == Dim::_3D) { + i1_prev(p) = i1(p); + dx1_prev(p) = dx1(p); + from_Xi_to_i_di(xp[0], i1(p), dx1(p)); + } - // update x2 & phi - if constexpr (D == Dim::_2D || D == Dim::_3D) { - i2_prev(p) = i2(p); - dx2_prev(p) = dx2(p); - from_Xi_to_i_di(xp[1], i2(p), dx2(p)); - if constexpr (D == Dim::_2D && M::PrtlDim == Dim::_3D) { - phi(p) = xp[2]; + // update x2 & phi + if constexpr (D == Dim::_2D || D == Dim::_3D) { + i2_prev(p) = i2(p); + dx2_prev(p) = dx2(p); + from_Xi_to_i_di(xp[1], i2(p), dx2(p)); + if constexpr (D == Dim::_2D && M::PrtlDim == Dim::_3D) { + phi(p) = xp[2]; + } } - } - // update x3 - if constexpr (D == Dim::_3D) { - i3_prev(p) = i3(p); - dx3_prev(p) = dx3(p); - from_Xi_to_i_di(xp[2], i3(p), dx3(p)); + // update x3 + if constexpr (D == Dim::_3D) { + i3_prev(p) = i3(p); + dx3_prev(p) = dx3(p); + from_Xi_to_i_di(xp[2], i3(p), dx3(p)); + } } boundaryConditions(p, xp); }