Merge pull request #3779 from lammps/fix-ttm-mod-arrays

Reorder fix ttm/mod 3d arrays to be consistent with fix ttm and fix ttm/grid
This commit is contained in:
Axel Kohlmeyer
2023-05-18 18:33:25 -04:00
committed by GitHub
4 changed files with 163 additions and 189 deletions

View File

@ -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
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: 1.781 ns/day, 13.479 hours/ns, 206.081 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
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

View File

@ -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
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: 5.393 ns/day, 4.450 hours/ns, 624.165 timesteps/s
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
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

View File

@ -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<Respa *>(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]);
else T_electron[iz][iy][ix] = T_electron_old[iz][iy][ix];
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;
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++];
}
/* ----------------------------------------------------------------------

View File

@ -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;