diff --git a/examples/ttm/log.20Apr22.ttm.mod.g++.1 b/examples/ttm/log.18May23.ttm.mod.g++.1 similarity index 70% rename from examples/ttm/log.20Apr22.ttm.mod.g++.1 rename to examples/ttm/log.18May23.ttm.mod.g++.1 index 6c0470617b2..b97e8ab0ea5 100644 --- a/examples/ttm/log.20Apr22.ttm.mod.g++.1 +++ b/examples/ttm/log.18May23.ttm.mod.g++.1 @@ -1,6 +1,4 @@ -LAMMPS (24 Mar 2022) -OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98) - using 1 OpenMP thread(s) per MPI task +LAMMPS (28 Mar 2023 - Development) units metal atom_style atomic boundary p p p @@ -15,7 +13,7 @@ mass 1 28.0855 create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 1 basis 5 1 basis 6 1 basis 7 1 basis 8 1 Created 8000 atoms using lattice units in orthogonal box = (0 0 0) to (54.309 54.309 54.309) - create_atoms CPU = 0.001 seconds + create_atoms CPU = 0.002 seconds pair_style sw pair_coeff * * Si.sw Si @@ -42,12 +40,12 @@ CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE Your simulation uses code contributions which should be cited: -- fix ttm/mod command: +- fix ttm/mod command: doi:10.1088/0953-8984/26/47/475401, doi:10.1002/ctpp.201310025 @article{Pisarev2014, author = {Pisarev, V. V. and Starikov, S. V.}, -title = {{Atomistic simulation of ion track formation in UO2.}}, -journal = {J.~Phys.:~Condens.~Matter}, +title = {Atomistic Simulation of Ion Track Formation in {UO$_2$}.}, +journal = {J.~Phys.\ Condens.\ Matter}, volume = {26}, number = {47}, pages = {475401}, @@ -56,8 +54,8 @@ year = {2014} @article{Norman2013, author = {Norman, G. E. and Starikov, S. V. and Stegailov, V. V. and Saitov, I. M. and Zhilyaev, P. A.}, -title = {{Atomistic Modeling of Warm Dense Matter in the Two-Temperature State}}, -journal = {Contrib.~Plasm.~Phys.}, +title = {Atomistic Modeling of Warm Dense Matter in the Two-Temperature State}, +journal = {Contrib.\ Plasma Phys.}, number = {2}, volume = {53}, pages = {129--139}, @@ -67,7 +65,7 @@ year = {2013} CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE Neighbor list info ... - update every 5 steps, delay 0 steps, check yes + update: every = 5 steps, delay = 0 steps, check = yes max neighbors/atom: 2000, page size: 100000 master list distance cutoff = 5.77118 ghost atom cutoff = 5.77118 @@ -81,30 +79,30 @@ Neighbor list info ... Per MPI rank memory allocation (min/avg/max) = 4.433 | 4.433 | 4.433 Mbytes Step Temp TotEng f_twotemp[1] f_twotemp[2] 0 0 -34692.79996100604 -52.79390940511979 0 - 100 2.004897156140836 -34690.27961013186 -55.34997305431884 0.01301140393178354 - 200 2.837118035232607 -34687.74741132015 -57.93445748841878 0.02696025968760173 - 300 4.263087164947482 -34684.98084093686 -60.75945453846786 0.02175636603841567 - 400 5.568003854939066 -34682.25271040963 -63.56896518300501 0.0300061848347275 - 500 6.225602451570786 -34679.49948952029 -66.40897551884576 0.02768827702656702 - 600 7.608847536264781 -34676.69728436362 -69.32060611557266 0.05579466731854093 - 700 9.049471241531297 -34674.00093206036 -72.10055094219446 0.004335980559879027 - 800 9.826796099683211 -34671.27720242751 -74.9501061086213 0.02371649678091513 - 900 11.8609224958918 -34668.35091308811 -77.98544170794551 0.004658649791374929 - 1000 13.88037467640968 -34665.35025858006 -81.16445160194114 0.07684078334464739 -Loop time of 4.85247 on 1 procs for 1000 steps with 8000 atoms - -Performance: 1.781 ns/day, 13.479 hours/ns, 206.081 timesteps/s -99.8% CPU use with 1 MPI tasks x 1 OpenMP threads + 100 2.004897156140836 -34690.27961013186 -55.3499730543189 0.01301140393178352 + 200 2.837118035232607 -34687.74741132015 -57.93445748841876 0.02696025968760173 + 300 4.263087164947482 -34684.98084093686 -60.75945453846793 0.02175636603841567 + 400 5.568003854939066 -34682.25271040963 -63.56896518300499 0.03000618483472749 + 500 6.225602451570786 -34679.49948952029 -66.40897551884574 0.02768827702656703 + 600 7.608847536264781 -34676.69728436362 -69.32060611557282 0.05579466731854091 + 700 9.049471241531297 -34674.00093206036 -72.10055094219462 0.004335980559879032 + 800 9.826796099683211 -34671.27720242751 -74.95010610862134 0.02371649678091515 + 900 11.8609224958918 -34668.35091308811 -77.98544170794545 0.004658649791374908 + 1000 13.88037467640968 -34665.35025858006 -81.16445160194111 0.07684078334464743 +Loop time of 2.48942 on 1 procs for 1000 steps with 8000 atoms + +Performance: 3.471 ns/day, 6.915 hours/ns, 401.700 timesteps/s, 3.214 Matom-step/s +100.0% CPU use with 1 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 4.1286 | 4.1286 | 4.1286 | 0.0 | 85.08 +Pair | 2.126 | 2.126 | 2.126 | 0.0 | 85.40 Neigh | 0 | 0 | 0 | 0.0 | 0.00 -Comm | 0.030972 | 0.030972 | 0.030972 | 0.0 | 0.64 -Output | 0.0026351 | 0.0026351 | 0.0026351 | 0.0 | 0.05 -Modify | 0.67848 | 0.67848 | 0.67848 | 0.0 | 13.98 -Other | | 0.01182 | | | 0.24 +Comm | 0.016147 | 0.016147 | 0.016147 | 0.0 | 0.65 +Output | 0.0013116 | 0.0013116 | 0.0013116 | 0.0 | 0.05 +Modify | 0.33864 | 0.33864 | 0.33864 | 0.0 | 13.60 +Other | | 0.007318 | | | 0.29 Nlocal: 8000 ave 8000 max 8000 min Histogram: 1 0 0 0 0 0 0 0 0 0 @@ -119,4 +117,4 @@ Total # of neighbors = 272000 Ave neighs/atom = 34 Neighbor list builds = 0 Dangerous builds = 0 -Total wall time: 0:00:04 +Total wall time: 0:00:02 diff --git a/examples/ttm/log.20Apr22.ttm.mod.g++.4 b/examples/ttm/log.18May23.ttm.mod.g++.4 similarity index 71% rename from examples/ttm/log.20Apr22.ttm.mod.g++.4 rename to examples/ttm/log.18May23.ttm.mod.g++.4 index fdf9b0cfb51..ea675c85941 100644 --- a/examples/ttm/log.20Apr22.ttm.mod.g++.4 +++ b/examples/ttm/log.18May23.ttm.mod.g++.4 @@ -1,6 +1,5 @@ -LAMMPS (24 Mar 2022) -OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98) - using 1 OpenMP thread(s) per MPI task +LAMMPS (28 Mar 2023 - Development) +WARNING: Using I/O redirection is unreliable with parallel runs. Better to use the -in switch to read input files. (../lammps.cpp:531) units metal atom_style atomic boundary p p p @@ -15,7 +14,7 @@ mass 1 28.0855 create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 1 basis 5 1 basis 6 1 basis 7 1 basis 8 1 Created 8000 atoms using lattice units in orthogonal box = (0 0 0) to (54.309 54.309 54.309) - create_atoms CPU = 0.000 seconds + create_atoms CPU = 0.001 seconds pair_style sw pair_coeff * * Si.sw Si @@ -42,12 +41,12 @@ CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE Your simulation uses code contributions which should be cited: -- fix ttm/mod command: +- fix ttm/mod command: doi:10.1088/0953-8984/26/47/475401, doi:10.1002/ctpp.201310025 @article{Pisarev2014, author = {Pisarev, V. V. and Starikov, S. V.}, -title = {{Atomistic simulation of ion track formation in UO2.}}, -journal = {J.~Phys.:~Condens.~Matter}, +title = {Atomistic Simulation of Ion Track Formation in {UO$_2$}.}, +journal = {J.~Phys.\ Condens.\ Matter}, volume = {26}, number = {47}, pages = {475401}, @@ -56,8 +55,8 @@ year = {2014} @article{Norman2013, author = {Norman, G. E. and Starikov, S. V. and Stegailov, V. V. and Saitov, I. M. and Zhilyaev, P. A.}, -title = {{Atomistic Modeling of Warm Dense Matter in the Two-Temperature State}}, -journal = {Contrib.~Plasm.~Phys.}, +title = {Atomistic Modeling of Warm Dense Matter in the Two-Temperature State}, +journal = {Contrib.\ Plasma Phys.}, number = {2}, volume = {53}, pages = {129--139}, @@ -67,7 +66,7 @@ year = {2013} CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE Neighbor list info ... - update every 5 steps, delay 0 steps, check yes + update: every = 5 steps, delay = 0 steps, check = yes max neighbors/atom: 2000, page size: 100000 master list distance cutoff = 5.77118 ghost atom cutoff = 5.77118 @@ -81,30 +80,30 @@ Neighbor list info ... Per MPI rank memory allocation (min/avg/max) = 3.436 | 3.436 | 3.436 Mbytes Step Temp TotEng f_twotemp[1] f_twotemp[2] 0 0 -34692.79996100361 -52.79390940511979 0 - 100 1.852689977101411 -34690.49204900486 -55.14271612882062 0.027261886765771 - 200 2.735750477179192 -34688.11139028054 -57.57110998717796 0.03387986355513582 - 300 3.931848271449558 -34685.54667417785 -60.18684521127226 0.02261256315262404 - 400 5.462009198576365 -34682.74455105668 -63.05420336037231 0.002402241637719583 - 500 6.267811692893873 -34679.96493887379 -65.93304222280049 0.02448378880218699 - 600 7.21148216150661 -34677.41455784726 -68.58391420045932 0.04114045759945373 - 700 8.84660534187052 -34674.40610468235 -71.68798344296847 0.0237298402743454 - 800 10.1748456457686 -34671.08749605772 -75.11943618276236 0.007538225788030298 - 900 11.27479036162859 -34668.4118066423 -77.92921692176769 0.02537529314475071 - 1000 13.26881394868076 -34665.56617589539 -80.91544540266329 0.03112665440209921 -Loop time of 1.60214 on 4 procs for 1000 steps with 8000 atoms - -Performance: 5.393 ns/day, 4.450 hours/ns, 624.165 timesteps/s -99.7% CPU use with 4 MPI tasks x 1 OpenMP threads + 100 1.852689977101411 -34690.49204900486 -55.14271612882064 0.02726188676577098 + 200 2.735750477179192 -34688.11139028054 -57.57110998717798 0.03387986355513584 + 300 3.931848271449558 -34685.54667417785 -60.18684521127231 0.02261256315262403 + 400 5.462009198576365 -34682.74455105668 -63.05420336037233 0.002402241637719578 + 500 6.267811692893873 -34679.96493887379 -65.93304222280051 0.02448378880218699 + 600 7.21148216150661 -34677.41455784726 -68.58391420045926 0.04114045759945374 + 700 8.84660534187052 -34674.40610468235 -71.68798344296859 0.02372984027434538 + 800 10.1748456457686 -34671.08749605772 -75.11943618276216 0.007538225788030307 + 900 11.27479036162859 -34668.4118066423 -77.92921692176756 0.02537529314475071 + 1000 13.26881394868076 -34665.56617589539 -80.91544540266317 0.03112665440209921 +Loop time of 0.995347 on 4 procs for 1000 steps with 8000 atoms + +Performance: 8.680 ns/day, 2.765 hours/ns, 1004.675 timesteps/s, 8.037 Matom-step/s +97.9% CPU use with 4 MPI tasks x no OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- -Pair | 1.0424 | 1.0558 | 1.0696 | 1.0 | 65.90 +Pair | 0.65351 | 0.6616 | 0.66783 | 0.8 | 66.47 Neigh | 0 | 0 | 0 | 0.0 | 0.00 -Comm | 0.05072 | 0.063773 | 0.079458 | 4.9 | 3.98 -Output | 0.0024362 | 0.0024703 | 0.0025297 | 0.1 | 0.15 -Modify | 0.47018 | 0.47332 | 0.48004 | 0.6 | 29.54 -Other | | 0.006786 | | | 0.42 +Comm | 0.041606 | 0.048314 | 0.056589 | 2.9 | 4.85 +Output | 0.0014609 | 0.0014742 | 0.0014968 | 0.0 | 0.15 +Modify | 0.27934 | 0.28016 | 0.28089 | 0.1 | 28.15 +Other | | 0.003798 | | | 0.38 Nlocal: 2000 ave 2000 max 2000 min Histogram: 4 0 0 0 0 0 0 0 0 0 diff --git a/src/EXTRA-FIX/fix_ttm_mod.cpp b/src/EXTRA-FIX/fix_ttm_mod.cpp index 1aa07b6781d..79af414f0a5 100644 --- a/src/EXTRA-FIX/fix_ttm_mod.cpp +++ b/src/EXTRA-FIX/fix_ttm_mod.cpp @@ -75,10 +75,8 @@ static constexpr double SHIFT = 0.0; FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - random(nullptr), nsum(nullptr), nsum_all(nullptr), - gfactor1(nullptr), gfactor2(nullptr), ratio(nullptr), flangevin(nullptr), - T_electron(nullptr), T_electron_old(nullptr), sum_vsq(nullptr), sum_mass_vsq(nullptr), - sum_vsq_all(nullptr), sum_mass_vsq_all(nullptr), net_energy_transfer(nullptr), + random(nullptr), gfactor1(nullptr), gfactor2(nullptr), ratio(nullptr), flangevin(nullptr), + T_electron(nullptr), T_electron_old(nullptr), net_energy_transfer(nullptr), net_energy_transfer_all(nullptr) { if (lmp->citeme) lmp->citeme->add(cite_fix_ttm_mod); @@ -167,20 +165,14 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) : // allocate 3d grid variables - memory->create(nsum,nxgrid,nygrid,nzgrid,"ttm/mod:nsum"); - memory->create(nsum_all,nxgrid,nygrid,nzgrid,"ttm/mod:nsum_all"); - memory->create(sum_vsq,nxgrid,nygrid,nzgrid,"ttm/mod:sum_vsq"); - memory->create(sum_mass_vsq,nxgrid,nygrid,nzgrid,"ttm/mod:sum_mass_vsq"); - memory->create(sum_vsq_all,nxgrid,nygrid,nzgrid,"ttm/mod:sum_vsq_all"); - memory->create(sum_mass_vsq_all,nxgrid,nygrid,nzgrid, - "ttm/mod:sum_mass_vsq_all"); - memory->create(T_electron_old,nxgrid,nygrid,nzgrid,"ttm/mod:T_electron_old"); - memory->create(T_electron_first,nxgrid,nygrid,nzgrid,"ttm/mod:T_electron_first"); - memory->create(T_electron,nxgrid,nygrid,nzgrid,"ttm/mod:T_electron"); - memory->create(net_energy_transfer,nxgrid,nygrid,nzgrid, + memory->create(T_electron_old,nzgrid,nygrid,nxgrid,"ttm/mod:T_electron_old"); + memory->create(T_electron_first,nzgrid,nygrid,nxgrid,"ttm/mod:T_electron_first"); + memory->create(T_electron,nzgrid,nygrid,nxgrid,"ttm/mod:T_electron"); + memory->create(net_energy_transfer,nzgrid,nygrid,nxgrid, "ttm/mod:net_energy_transfer"); - memory->create(net_energy_transfer_all,nxgrid,nygrid,nzgrid, + memory->create(net_energy_transfer_all,nzgrid,nygrid,nxgrid, "ttm/mod:net_energy_transfer_all"); + flangevin = nullptr; grow_arrays(atom->nmax); @@ -203,10 +195,10 @@ FixTTMMod::FixTTMMod(LAMMPS *lmp, int narg, char **arg) : // initialize electron temperatures on grid int ix,iy,iz; - for (ix = 0; ix < nxgrid; ix++) + for (iz = 0; iz < nzgrid; iz++) for (iy = 0; iy < nygrid; iy++) - for (iz = 0; iz < nzgrid; iz++) - T_electron[ix][iy][iz] = tinit; + for (ix = 0; ix < nxgrid; ix++) + T_electron[iz][iy][ix] = tinit; // if specified, read initial electron temperatures from file @@ -221,18 +213,13 @@ FixTTMMod::~FixTTMMod() delete[] gfactor1; delete[] gfactor2; - memory->destroy(nsum); - memory->destroy(nsum_all); - memory->destroy(sum_vsq); - memory->destroy(sum_mass_vsq); - memory->destroy(sum_vsq_all); - memory->destroy(sum_mass_vsq_all); memory->destroy(T_electron_first); memory->destroy(T_electron_old); memory->destroy(T_electron); - memory->destroy(flangevin); memory->destroy(net_energy_transfer); memory->destroy(net_energy_transfer_all); + + memory->destroy(flangevin); } /* ---------------------------------------------------------------------- */ @@ -265,10 +252,10 @@ void FixTTMMod::init() sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v; } - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) - net_energy_transfer_all[ix][iy][iz] = 0; + for (int ix = 0; ix < nxgrid; ix++) + net_energy_transfer_all[iz][iy][ix] = 0; if (utils::strmatch(update->integrate_style,"^respa")) nlevels_respa = (dynamic_cast(update->integrate))->nlevels; @@ -320,10 +307,10 @@ void FixTTMMod::post_force(int /*vflag*/) while (iy < 0) iy += nygrid; while (iz < 0) iz += nzgrid; - if (T_electron[ix][iy][iz] < 0) + if (T_electron[iz][iy][ix] < 0) error->all(FLERR,"Electronic temperature dropped below zero"); - double tsqrt = sqrt(T_electron[ix][iy][iz]); + double tsqrt = sqrt(T_electron[iz][iy][ix]); gamma1 = gfactor1[type[i]]; double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; @@ -342,20 +329,14 @@ void FixTTMMod::post_force(int /*vflag*/) if (right_x == nxgrid) right_x = 0; if (right_y == nygrid) right_y = 0; if (right_z == nzgrid) right_z = 0; - int left_x = ix - 1; - int left_y = iy - 1; - int left_z = iz - 1; - if (left_x == -1) left_x = nxgrid - 1; - if (left_y == -1) left_y = nygrid - 1; - if (left_z == -1) left_z = nzgrid - 1; - double T_i = T_electron[ix][iy][iz]; - double T_ir = T_electron[right_x][iy][iz]; - double T_iu = T_electron[ix][right_y][iz]; - double T_if = T_electron[ix][iy][right_z]; - double C_i = el_properties(T_electron[ix][iy][iz]).el_heat_capacity; - double C_ir = el_properties(T_electron[right_x][iy][iz]).el_heat_capacity; - double C_iu = el_properties(T_electron[ix][right_y][iz]).el_heat_capacity; - double C_if = el_properties(T_electron[ix][iy][right_z]).el_heat_capacity; + double T_i = T_electron[iz][iy][ix]; + double T_ir = T_electron[iz][iy][right_x]; + double T_iu = T_electron[iz][right_y][ix]; + double T_if = T_electron[right_z][iy][ix]; + double C_i = el_properties(T_electron[iz][iy][ix]).el_heat_capacity; + double C_ir = el_properties(T_electron[iz][iy][right_x]).el_heat_capacity; + double C_iu = el_properties(T_electron[iz][right_y][ix]).el_heat_capacity; + double C_if = el_properties(T_electron[right_z][iy][ix]).el_heat_capacity; double diff_x = (x_at - x_surf)*(x_at - x_surf); diff_x = pow(diff_x,0.5); double len_factor = diff_x/(diff_x+free_path); @@ -587,7 +568,7 @@ void FixTTMMod::read_electron_temperatures(const std::string &filename) if (comm->me == 0) { int ***T_initial_set; - memory->create(T_initial_set,nxgrid,nygrid,nzgrid,"ttm/mod:T_initial_set"); + memory->create(T_initial_set,nzgrid,nygrid,nxgrid,"ttm/mod:T_initial_set"); memset(&T_initial_set[0][0][0],0,ngridtotal*sizeof(int)); // read initial electron temperature values from file @@ -614,8 +595,8 @@ void FixTTMMod::read_electron_temperatures(const std::string &filename) if (T_tmp < 0.0) throw TokenizerException("Fix ttm electron temperatures must be > 0.0",""); - T_electron[ix][iy][iz] = T_tmp; - T_initial_set[ix][iy][iz] = 1; + T_electron[iz][iy][ix] = T_tmp; + T_initial_set[iz][iy][ix] = 1; } } catch (std::exception &e) { error->one(FLERR, e.what()); @@ -626,7 +607,7 @@ void FixTTMMod::read_electron_temperatures(const std::string &filename) for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) for (int ix = 0; ix < nxgrid; ix++) - if (T_initial_set[ix][iy][iz] == 0) + if (T_initial_set[iz][iy][ix] == 0) error->all(FLERR,"Fix ttm/mod infile did not set all temperatures"); memory->destroy(T_initial_set); @@ -656,9 +637,9 @@ void FixTTMMod::write_electron_temperatures(const std::string &filename) for (ix = 0; ix < nxgrid; ix++) for (iy = 0; iy < nygrid; iy++) for (iz = 0; iz < nzgrid; iz++) { - if (movsur == 1 && T_electron[ix][iy][iz] == 0.0) - T_electron[ix][iy][iz] = electron_temperature_min; - fprintf(fp,"%d %d %d %20.16g\n",ix+1,iy+1,iz+1,T_electron[ix][iy][iz]); + if (movsur == 1 && T_electron[iz][iy][ix] == 0.0) + T_electron[iz][iy][ix] = electron_temperature_min; + fprintf(fp,"%d %d %d %20.16g\n",ix+1,iy+1,iz+1,T_electron[iz][iy][ix]); } fclose(fp); @@ -678,6 +659,9 @@ el_heat_capacity_thermal_conductivity FixTTMMod::el_properties(double T_e) properties.el_thermal_conductivity = el_th_diff*properties.el_heat_capacity; // thermal conductivity return properties; } + +/* ---------------------------------------------------------------------- */ + double FixTTMMod::el_sp_heat_integral(double T_e) { double T_temp = T_e/1000.0, T_reduced = T_damp*T_temp; @@ -700,20 +684,21 @@ void FixTTMMod::end_of_step() int nlocal = atom->nlocal; if (movsur == 1) { - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) { - double TTT = T_electron[ix][iy][iz]; + for (int ix = 0; ix < nxgrid; ix++) { + double TTT = T_electron[iz][iy][ix]; if (TTT > 0) { if (ix < t_surface_l) t_surface_l = ix; } } } - for (int ix = 0; ix < nxgrid; ix++) + + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) - net_energy_transfer[ix][iy][iz] = 0; + for (int ix = 0; ix < nxgrid; ix++) + net_energy_transfer[iz][iy][ix] = 0; for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -731,7 +716,7 @@ void FixTTMMod::end_of_step() while (iz < 0) iz += nzgrid; if (ix >= t_surface_l) { if (ix < surface_r) - net_energy_transfer[ix][iy][iz] += + net_energy_transfer[iz][iy][ix] += (flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] + flangevin[i][2]*v[i][2]); } @@ -747,18 +732,16 @@ void FixTTMMod::end_of_step() double del_vol = dx*dy*dz; double el_specific_heat = 0.0; double el_thermal_conductivity = el_properties(electron_temperature_min).el_thermal_conductivity; - for (int ix = 0; ix < nxgrid; ix++) + + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) - { - if (el_properties(T_electron[ix][iy][iz]).el_thermal_conductivity > el_thermal_conductivity) - el_thermal_conductivity = el_properties(T_electron[ix][iy][iz]).el_thermal_conductivity; - if (el_specific_heat > 0.0) - { - if ((T_electron[ix][iy][iz] > 0.0) && (el_properties(T_electron[ix][iy][iz]).el_heat_capacity < el_specific_heat)) - el_specific_heat = el_properties(T_electron[ix][iy][iz]).el_heat_capacity; - } - else if (T_electron[ix][iy][iz] > 0.0) el_specific_heat = el_properties(T_electron[ix][iy][iz]).el_heat_capacity; + for (int ix = 0; ix < nxgrid; ix++) { + if (el_properties(T_electron[iz][iy][ix]).el_thermal_conductivity > el_thermal_conductivity) + el_thermal_conductivity = el_properties(T_electron[iz][iy][ix]).el_thermal_conductivity; + if (el_specific_heat > 0.0) { + if ((T_electron[iz][iy][ix] > 0.0) && (el_properties(T_electron[iz][iy][ix]).el_heat_capacity < el_specific_heat)) + el_specific_heat = el_properties(T_electron[iz][iy][ix]).el_heat_capacity; + } else if (T_electron[iz][iy][ix] > 0.0) el_specific_heat = el_properties(T_electron[iz][iy][ix]).el_heat_capacity; } // num_inner_timesteps = # of inner steps (thermal solves) // required this MD step to maintain a stable explicit solve @@ -767,17 +750,16 @@ void FixTTMMod::end_of_step() double inner_dt = update->dt; double stability_criterion = 0.0; - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) - T_electron_first[ix][iy][iz] = - T_electron[ix][iy][iz]; + for (int ix = 0; ix < nxgrid; ix++) + T_electron_first[iz][iy][ix] = T_electron[iz][iy][ix]; + do { - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) - T_electron[ix][iy][iz] = - T_electron_first[ix][iy][iz]; + for (int ix = 0; ix < nxgrid; ix++) + T_electron[iz][iy][ix] = T_electron_first[iz][iy][ix]; stability_criterion = 1.0 - 2.0*inner_dt/el_specific_heat * @@ -792,16 +774,15 @@ void FixTTMMod::end_of_step() error->warning(FLERR,"Too many inner timesteps in fix ttm/mod"); for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps; ith_inner_timestep++) { - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) - T_electron_old[ix][iy][iz] = - T_electron[ix][iy][iz]; + for (int ix = 0; ix < nxgrid; ix++) + T_electron_old[iz][iy][ix] = T_electron[iz][iy][ix]; // compute new electron T profile duration = duration + inner_dt; - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) { + for (int ix = 0; ix < nxgrid; ix++) { int right_x = ix + 1; int right_y = iy + 1; int right_z = iz + 1; @@ -821,50 +802,50 @@ void FixTTMMod::end_of_step() if (duration < width) { if (ix >= t_surface_l) mult_factor = (intensity/(dx*skin_layer_d))*exp((-1.0)*(ix_d - surface_d)/skin_layer_d); } - if (ix < t_surface_l) net_energy_transfer_all[ix][iy][iz] = 0.0; + if (ix < t_surface_l) net_energy_transfer_all[iz][iy][ix] = 0.0; double cr_vac = 1; - if (T_electron_old[ix][iy][iz] == 0) cr_vac = 0; + if (T_electron_old[iz][iy][ix] == 0) cr_vac = 0; double cr_v_l_x = 1; - if (T_electron_old[left_x][iy][iz] == 0) cr_v_l_x = 0; + if (T_electron_old[iz][iy][left_x] == 0) cr_v_l_x = 0; double cr_v_r_x = 1; - if (T_electron_old[right_x][iy][iz] == 0) cr_v_r_x = 0; + if (T_electron_old[iz][iy][right_x] == 0) cr_v_r_x = 0; double cr_v_l_y = 1; - if (T_electron_old[ix][left_y][iz] == 0) cr_v_l_y = 0; + if (T_electron_old[iz][left_y][ix] == 0) cr_v_l_y = 0; double cr_v_r_y = 1; - if (T_electron_old[ix][right_y][iz] == 0) cr_v_r_y = 0; + if (T_electron_old[iz][right_y][ix] == 0) cr_v_r_y = 0; double cr_v_l_z = 1; - if (T_electron_old[ix][iy][left_z] == 0) cr_v_l_z = 0; + if (T_electron_old[left_z][iy][ix] == 0) cr_v_l_z = 0; double cr_v_r_z = 1; - if (T_electron_old[ix][iy][right_z] == 0) cr_v_r_z = 0; + if (T_electron_old[right_z][iy][ix] == 0) cr_v_r_z = 0; if (cr_vac != 0) { - T_electron[ix][iy][iz] = - T_electron_old[ix][iy][iz] + - inner_dt/el_properties(T_electron_old[ix][iy][iz]).el_heat_capacity * - ((cr_v_r_x*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[right_x][iy][iz]/2.0).el_thermal_conductivity* - (T_electron_old[right_x][iy][iz]-T_electron_old[ix][iy][iz])/dx - - cr_v_l_x*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[left_x][iy][iz]/2.0).el_thermal_conductivity* - (T_electron_old[ix][iy][iz]-T_electron_old[left_x][iy][iz])/dx)/dx + - (cr_v_r_y*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[ix][right_y][iz]/2.0).el_thermal_conductivity* - (T_electron_old[ix][right_y][iz]-T_electron_old[ix][iy][iz])/dy - - cr_v_l_y*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[ix][left_y][iz]/2.0).el_thermal_conductivity* - (T_electron_old[ix][iy][iz]-T_electron_old[ix][left_y][iz])/dy)/dy + - (cr_v_r_z*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[ix][iy][right_z]/2.0).el_thermal_conductivity* - (T_electron_old[ix][iy][right_z]-T_electron_old[ix][iy][iz])/dz - - cr_v_l_z*el_properties(T_electron_old[ix][iy][iz]/2.0+T_electron_old[ix][iy][left_z]/2.0).el_thermal_conductivity* - (T_electron_old[ix][iy][iz]-T_electron_old[ix][iy][left_z])/dz)/dz); - T_electron[ix][iy][iz]+=inner_dt/el_properties(T_electron[ix][iy][iz]).el_heat_capacity* + T_electron[iz][iy][ix] = + T_electron_old[iz][iy][ix] + + inner_dt/el_properties(T_electron_old[iz][iy][ix]).el_heat_capacity * + ((cr_v_r_x*el_properties(T_electron_old[iz][iy][ix]/2.0+T_electron_old[iz][iy][right_x]/2.0).el_thermal_conductivity* + (T_electron_old[iz][iy][right_x]-T_electron_old[iz][iy][ix])/dx - + cr_v_l_x*el_properties(T_electron_old[iz][iy][ix]/2.0+T_electron_old[iz][iy][left_x]/2.0).el_thermal_conductivity* + (T_electron_old[iz][iy][ix]-T_electron_old[iz][iy][left_x])/dx)/dx + + (cr_v_r_y*el_properties(T_electron_old[iz][iy][ix]/2.0+T_electron_old[iz][right_y][ix]/2.0).el_thermal_conductivity* + (T_electron_old[iz][right_y][ix]-T_electron_old[iz][iy][ix])/dy - + cr_v_l_y*el_properties(T_electron_old[iz][iy][ix]/2.0+T_electron_old[iz][left_y][ix]/2.0).el_thermal_conductivity* + (T_electron_old[iz][iy][ix]-T_electron_old[iz][left_y][ix])/dy)/dy + + (cr_v_r_z*el_properties(T_electron_old[iz][iy][ix]/2.0+T_electron_old[right_z][iy][ix]/2.0).el_thermal_conductivity* + (T_electron_old[right_z][iy][ix]-T_electron_old[iz][iy][ix])/dz - + cr_v_l_z*el_properties(T_electron_old[iz][iy][ix]/2.0+T_electron_old[left_z][iy][ix]/2.0).el_thermal_conductivity* + (T_electron_old[iz][iy][ix]-T_electron_old[left_z][iy][ix])/dz)/dz); + T_electron[iz][iy][ix]+=inner_dt/el_properties(T_electron[iz][iy][ix]).el_heat_capacity* (mult_factor - - net_energy_transfer_all[ix][iy][iz]/del_vol); + net_energy_transfer_all[iz][iy][ix]/del_vol); } - else T_electron[ix][iy][iz] = - T_electron_old[ix][iy][iz]; - if ((T_electron[ix][iy][iz] > 0.0) && (T_electron[ix][iy][iz] < electron_temperature_min)) - T_electron[ix][iy][iz] = T_electron[ix][iy][iz] + 0.5*(electron_temperature_min - T_electron[ix][iy][iz]); - - if (el_properties(T_electron[ix][iy][iz]).el_thermal_conductivity > el_thermal_conductivity) - el_thermal_conductivity = el_properties(T_electron[ix][iy][iz]).el_thermal_conductivity; - if ((T_electron[ix][iy][iz] > 0.0) && (el_properties(T_electron[ix][iy][iz]).el_heat_capacity < el_specific_heat)) - el_specific_heat = el_properties(T_electron[ix][iy][iz]).el_heat_capacity; + else T_electron[iz][iy][ix] = T_electron_old[iz][iy][ix]; + + if ((T_electron[iz][iy][ix] > 0.0) && (T_electron[iz][iy][ix] < electron_temperature_min)) + T_electron[iz][iy][ix] = T_electron[iz][iy][ix] + 0.5*(electron_temperature_min - T_electron[iz][iy][ix]); + + if (el_properties(T_electron[iz][iy][ix]).el_thermal_conductivity > el_thermal_conductivity) + el_thermal_conductivity = el_properties(T_electron[iz][iy][ix]).el_thermal_conductivity; + if ((T_electron[iz][iy][ix] > 0.0) && (el_properties(T_electron[iz][iy][ix]).el_heat_capacity < el_specific_heat)) + el_specific_heat = el_properties(T_electron[iz][iy][ix]).el_heat_capacity; } } stability_criterion = 1.0 - @@ -913,12 +894,11 @@ double FixTTMMod::compute_vector(int n) double dz = domain->zprd/nzgrid; double del_vol = dx*dy*dz; - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) { - e_energy += el_sp_heat_integral(T_electron[ix][iy][iz])*del_vol; - transfer_energy += - net_energy_transfer_all[ix][iy][iz]*update->dt; + for (int ix = 0; ix < nxgrid; ix++) { + e_energy += el_sp_heat_integral(T_electron[iz][iy][ix])*del_vol; + transfer_energy += net_energy_transfer_all[iz][iy][ix]*update->dt; } if (n == 0) return e_energy; @@ -933,7 +913,7 @@ double FixTTMMod::compute_vector(int n) void FixTTMMod::write_restart(FILE *fp) { double *rlist; - memory->create(rlist,nxgrid*nygrid*nzgrid+4,"ttm/mod:rlist"); + memory->create(rlist,nzgrid*nygrid*nxgrid + 4,"ttm/mod:rlist"); int n = 0; rlist[n++] = nxgrid; @@ -941,10 +921,10 @@ void FixTTMMod::write_restart(FILE *fp) rlist[n++] = nzgrid; rlist[n++] = seed; - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) - rlist[n++] = T_electron[ix][iy][iz]; + for (int ix = 0; ix < nxgrid; ix++) + rlist[n++] = T_electron[iz][iy][ix]; if (comm->me == 0) { int size = n * sizeof(double); @@ -982,10 +962,10 @@ void FixTTMMod::restart(char *buf) // restore global frid values - for (int ix = 0; ix < nxgrid; ix++) + for (int iz = 0; iz < nzgrid; iz++) for (int iy = 0; iy < nygrid; iy++) - for (int iz = 0; iz < nzgrid; iz++) - T_electron[ix][iy][iz] = rlist[n++]; + for (int ix = 0; ix < nxgrid; ix++) + T_electron[iz][iy][ix] = rlist[n++]; } /* ---------------------------------------------------------------------- diff --git a/src/EXTRA-FIX/fix_ttm_mod.h b/src/EXTRA-FIX/fix_ttm_mod.h index 28c32096976..469a641d567 100644 --- a/src/EXTRA-FIX/fix_ttm_mod.h +++ b/src/EXTRA-FIX/fix_ttm_mod.h @@ -64,11 +64,8 @@ class FixTTMMod : public Fix { int nxgrid, nygrid, nzgrid; int ngridtotal; - int ***nsum, ***nsum_all; double *gfactor1, *gfactor2, *ratio, **flangevin; double ***T_electron, ***T_electron_old, ***T_electron_first; - double ***sum_vsq, ***sum_mass_vsq; - double ***sum_vsq_all, ***sum_mass_vsq_all; double ***net_energy_transfer, ***net_energy_transfer_all; double gamma_p, gamma_s, v_0, v_0_sq;