Fixing errors in add heat

This commit is contained in:
jtclemm
2024-08-09 15:32:52 -06:00
parent 07ac930733
commit a7c6b0ee42
2 changed files with 13 additions and 14 deletions

View File

@ -18,11 +18,11 @@ Syntax
*constant* args = *rate* *constant* args = *rate*
*rate* = rate of heat flow (energy/time units) *rate* = rate of heat flow (energy/time units)
*linear* args = :math:`T_{target}` *k* *linear* args = *Ttarget* *k*
:math:`T_{target}` = target temperature (temperature units) *Ttarget* = target temperature (temperature units)
*k* = prefactor (energy/(time*temperature) units) *k* = prefactor (energy/(time*temperature) units)
*quartic* args = :math:`T_{target}` *k* *quartic* args = *Ttarget* *k*
:math:`T_{target}` = target temperature (temperature units) *Ttarget* = target temperature (temperature units)
*k* = prefactor (energy/(time*temperature^4) units) *k* = prefactor (energy/(time*temperature^4) units)
* zero or more keyword/value pairs may be appended to args * zero or more keyword/value pairs may be appended to args
@ -45,9 +45,9 @@ Examples
Description Description
""""""""""" """""""""""
This fix adds heat to particles with the temperature attribute every timestep. This fix adds heat to particles with the temperature attribute every timestep
Note that this is an internal temperature of a particle intended for use with at a given rate. Note that this is an internal temperature of a particle intended
non-atomistic models like the discrete element method. for use with non-atomistic models like the discrete element method.
For the *constant* style, heat is added at the specified rate. For the *linear* style, For the *constant* style, heat is added at the specified rate. For the *linear* style,
heat is added at a rate of :math:`k (T_{target} - T)` where :math:`k` is the heat is added at a rate of :math:`k (T_{target} - T)` where :math:`k` is the
@ -96,7 +96,8 @@ only enabled if LAMMPS was built with that package.
See the :doc:`Build package <Build_package>` page for more info. See the :doc:`Build package <Build_package>` page for more info.
This fix requires that atoms store temperature and heat flow This fix requires that atoms store temperature and heat flow
as defined by the :doc:`fix property/atom <fix_property_atom>` command. as defined by the :doc:`fix property/atom <fix_property_atom>` command or
included in certain atom styles, such as atom_style rheo/thermal.
Related commands Related commands
"""""""""""""""" """"""""""""""""

View File

@ -21,7 +21,6 @@
#include "error.h" #include "error.h"
#include "input.h" #include "input.h"
#include "memory.h" #include "memory.h"
#include "update.h"
#include "variable.h" #include "variable.h"
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -125,7 +124,6 @@ void FixAddHeat::post_force(int /*vflag*/)
int *mask = atom->mask; int *mask = atom->mask;
double *heatflow = atom->heatflow; double *heatflow = atom->heatflow;
double *temperature = atom->temperature; double *temperature = atom->temperature;
double dtinv = 1.0 / update->dt;
if (vstyle == ATOM) { if (vstyle == ATOM) {
if (atom->nmax > maxatom) { if (atom->nmax > maxatom) {
@ -142,7 +140,7 @@ void FixAddHeat::post_force(int /*vflag*/)
if (mask[i] & groupbit) if (mask[i] & groupbit)
heatflow[i] = 0.0; heatflow[i] = 0.0;
double vtmp, dt; double vtmp;
if (vstyle == CONSTANT) vtmp = value; if (vstyle == CONSTANT) vtmp = value;
if (vstyle == EQUAL) vtmp = input->variable->compute_equal(var); if (vstyle == EQUAL) vtmp = input->variable->compute_equal(var);
for (int i = 0; i < atom->nlocal; i++) { for (int i = 0; i < atom->nlocal; i++) {
@ -150,11 +148,11 @@ void FixAddHeat::post_force(int /*vflag*/)
if (vstyle == ATOM) vtmp = vatom[i]; if (vstyle == ATOM) vtmp = vatom[i];
if (style == ADD) { if (style == ADD) {
heatflow[i] += dtinv * vtmp; heatflow[i] += vtmp;
} else if (style == LINEAR) { } else if (style == LINEAR) {
heatflow[i] += dtinv * prefactor * (vtmp - temperature[i]); heatflow[i] += prefactor * (vtmp - temperature[i]);
} else if (style == QUARTIC) { } else if (style == QUARTIC) {
heatflow[i] += dtinv * prefactor * (pow(vtmp, 4.0) - pow(temperature[i], 4.0)); heatflow[i] += prefactor * (pow(vtmp, 4.0) - pow(temperature[i], 4.0));
} }
} }
} }