/* ----------------------------------------------------------------------
 *
 *                    *** Smooth Mach Dynamics ***
 *
 * This file is part of the USER-SMD package for LAMMPS.
 * Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
 * Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
 * Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
 *
 * ----------------------------------------------------------------------- */

/* ----------------------------------------------------------------------
 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
 http://lammps.sandia.gov, Sandia National Laboratories
 Steve Plimpton, sjplimp@sandia.gov

 Copyright (2003) Sandia Corporation.  Under the terms of Contract
 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
 certain rights in this software.  This software is distributed under
 the GNU General Public License.

 See the README file in the top-level LAMMPS directory.
 ------------------------------------------------------------------------- */

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "fix_smd_adjust_dt.h"
#include "atom.h"
#include "update.h"
#include "integrate.h"
#include "domain.h"
#include "lattice.h"
#include "force.h"
#include "pair.h"
#include "modify.h"
#include "fix.h"
#include "output.h"
#include "dump.h"
#include "comm.h"
#include "error.h"

using namespace LAMMPS_NS;
using namespace FixConst;

#define BIG 1.0e20

/* ---------------------------------------------------------------------- */

FixSMDTlsphDtReset::FixSMDTlsphDtReset(LAMMPS *lmp, int narg, char **arg) :
		Fix(lmp, narg, arg) {
	if (narg != 4)
		error->all(FLERR, "Illegal fix smd/adjust_dt command");

	// set time_depend, else elapsed time accumulation can be messed up

	time_depend = 1;
	scalar_flag = 1;
	vector_flag = 1;
	size_vector = 2;
	global_freq = 1;
	extscalar = 0;
	extvector = 0;
	restart_global = 1; // this fix stores global (i.e., not per-atom) info: elaspsed time

	safety_factor = atof(arg[3]);

	// initializations
	t_elapsed = 0.0;
}

/* ---------------------------------------------------------------------- */

int FixSMDTlsphDtReset::setmask() {
	int mask = 0;
	mask |= INITIAL_INTEGRATE;
	mask |= END_OF_STEP;
	return mask;
}

/* ---------------------------------------------------------------------- */

void FixSMDTlsphDtReset::init() {
	dt = update->dt;
}

/* ---------------------------------------------------------------------- */

void FixSMDTlsphDtReset::setup(int vflag) {
	end_of_step();
}

/* ---------------------------------------------------------------------- */

void FixSMDTlsphDtReset::initial_integrate(int vflag) {

	//printf("in adjust_dt: dt = %20.10f\n", update->dt);

	t_elapsed += update->dt;
}

/* ---------------------------------------------------------------------- */

void FixSMDTlsphDtReset::end_of_step() {
	double dtmin = BIG;
	int itmp = 0;

	/*
	 * extract minimum CFL timestep from TLSPH and ULSPH pair styles
	 */

	double *dtCFL_TLSPH = (double *) force->pair->extract("smd/tlsph/dtCFL_ptr", itmp);
	double *dtCFL_ULSPH = (double *) force->pair->extract("smd/ulsph/dtCFL_ptr", itmp);
	double *dt_TRI = (double *) force->pair->extract("smd/tri_surface/stable_time_increment_ptr", itmp);
	double *dt_HERTZ = (double *) force->pair->extract("smd/hertz/stable_time_increment_ptr", itmp);
	double *dt_PERI_IPMB = (double *) force->pair->extract("smd/peri_ipmb/stable_time_increment_ptr", itmp);

	if ((dtCFL_TLSPH == NULL) && (dtCFL_ULSPH == NULL) && (dt_TRI == NULL) && (dt_HERTZ == NULL)
			&& (dt_PERI_IPMB == NULL)) {
		error->all(FLERR, "fix smd/adjust_dt failed to access a valid dtCFL");
	}

	if (dtCFL_TLSPH != NULL) {
		dtmin = MIN(dtmin, *dtCFL_TLSPH);
	}

	if (dtCFL_ULSPH != NULL) {
		dtmin = MIN(dtmin, *dtCFL_ULSPH);
	}

	if (dt_TRI != NULL) {
		dtmin = MIN(dtmin, *dt_TRI);
	}

	if (dt_HERTZ != NULL) {
		dtmin = MIN(dtmin, *dt_HERTZ);
	}

	if (dt_PERI_IPMB != NULL) {
		dtmin = MIN(dtmin, *dt_PERI_IPMB);
	}

//	double **v = atom->v;
//	double **f = atom->f;
//	double *rmass = atom->rmass;
//	double *radius = atom->radius;
//	int *mask = atom->mask;
//	int nlocal = atom->nlocal;
//	double dtv, dtf, dtsq;
//	double vsq, fsq, massinv, xmax;
//	double delx, dely, delz, delr;

//	for (int i = 0; i < nlocal; i++) {
//		if (mask[i] & groupbit) {
//			xmax = 0.005 * radius[i];
//			massinv = 1.0 / rmass[i];
//			vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
//			fsq = f[i][0] * f[i][0] + f[i][1] * f[i][1] + f[i][2] * f[i][2];
//			dtv = dtf = BIG;
//			if (vsq > 0.0)
//				dtv = xmax / sqrt(vsq);
//			if (fsq > 0.0)
//				dtf = sqrt(2.0 * xmax / (sqrt(fsq) * massinv));
//			dt = MIN(dtv, dtf);
//			dtmin = MIN(dtmin, dt);
//			dtsq = dt * dt;
//			delx = dt * v[i][0] + 0.5 * dtsq * massinv * f[i][0];
//			dely = dt * v[i][1] + 0.5 * dtsq * massinv * f[i][1];
//			delz = dt * v[i][2] + 0.5 * dtsq * massinv * f[i][2];
//			delr = sqrt(delx * delx + dely * dely + delz * delz);
//			if (delr > xmax)
//				dt *= xmax / delr;
//			dtmin = MIN(dtmin, dt);
//
////			xmax = 0.05 * radius[i];
////			massinv = 1.0 / rmass[i];
////			fsq = f[i][0] * f[i][0] + f[i][1] * f[i][1] + f[i][2] * f[i][2];
////			dtf = BIG;
////			if (fsq > 0.0)
////				dtf = sqrt(2.0 * xmax / (sqrt(fsq) * massinv));
////			dtmin = MIN(dtmin, dtf);
//		}
//	}

	dtmin *= safety_factor; // apply safety factor
	MPI_Allreduce(&dtmin, &dt, 1, MPI_DOUBLE, MPI_MIN, world);

	if (update->ntimestep == 0) {
		dt = 1.0e-16;
	}

	//printf("dtmin is now: %f, dt is now%f\n", dtmin, dt);


	update->dt = dt;
	if (force->pair)
		force->pair->reset_dt();
	for (int i = 0; i < modify->nfix; i++)
		modify->fix[i]->reset_dt();
}

/* ---------------------------------------------------------------------- */

double FixSMDTlsphDtReset::compute_scalar() {
	return t_elapsed;
}

/* ----------------------------------------------------------------------
 pack entire state of Fix into one write
 ------------------------------------------------------------------------- */

void FixSMDTlsphDtReset::write_restart(FILE *fp) {
	int n = 0;
	double list[1];
	list[n++] = t_elapsed;

	if (comm->me == 0) {
		int size = n * sizeof(double);
		fwrite(&size, sizeof(int), 1, fp);
		fwrite(list, sizeof(double), n, fp);
	}
}

/* ----------------------------------------------------------------------
 use state info from restart file to restart the Fix
 ------------------------------------------------------------------------- */

void FixSMDTlsphDtReset::restart(char *buf) {
	int n = 0;
	double *list = (double *) buf;
	t_elapsed = list[n++];
}

