applications/hydrodynamics/Hydro.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
#include <iostream>
#include <fstream>
#include <string>

#include "applications/hydrodynamics/HydroCmd.h"
#include "applications/hydrodynamics/HydrodynamicsModel.hpp"
#include "applications/hydrodynamics/HydrodynamicsModelCreator.hpp"
#include "applications/hydrodynamics/HydrodynamicsModelFactory.hpp"
#include "applications/hydrodynamics/AnalyticalModel.hpp"
#include "applications/hydrodynamics/BeadModel.hpp"
#include "applications/hydrodynamics/RoughShell.hpp"
#include "applications/hydrodynamics/ShapeBuilder.hpp"
#include "brains/Register.hpp"
#include "brains/SimCreator.hpp"
#include "brains/SimInfo.hpp"
#include "utils/StringUtils.hpp"
#include "utils/simError.h"
#include "utils/MemoryUtils.hpp"
using namespace OpenMD;

struct SDShape{
    StuntDouble* sd;
    Shape* shape;
};
void registerHydrodynamicsModels();
void writeHydroProps(std::ostream& os);

int main(int argc, char* argv[]){
  registerHydrodynamicsModels();
  
  gengetopt_args_info args_info;
  std::string dumpFileName;
  std::string mdFileName;
  std::string prefix;
  
  //parse the command line option
  if (cmdline_parser (argc, argv, &args_info) != 0) {
    exit(1) ;
  }
  
  //get the dumpfile name and meta-data file name
  if (args_info.input_given){
    dumpFileName = args_info.input_arg;
  } else {
    std::cerr << "Does not have input file name" << std::endl;
    exit(1);
  }
  
  if (args_info.output_given){
    prefix = args_info.output_arg;
  } else {
    prefix = getPrefix(dumpFileName);    
  }
  std::string outputFilename = prefix + ".diff";
    
  //parse md file and set up the system
  SimCreator creator;
  SimInfo* info = creator.createSim(dumpFileName, true);
    
  SimInfo::MoleculeIterator mi;
  Molecule* mol;
  Molecule::IntegrableObjectIterator  ii;
  StuntDouble* integrableObject;
  Mat3x3d identMat;
  identMat(0,0) = 1.0;
  identMat(1,1) = 1.0;
  identMat(2,2) = 1.0;

  Globals* simParams = info->getSimParams();
  RealType temperature;
  RealType viscosity;

  if (simParams->haveViscosity()) {
    viscosity = simParams->getViscosity();
  } else {
    sprintf(painCave.errMsg, "viscosity must be set\n");
    painCave.isFatal = 1;
    simError();  
  }

  if (simParams->haveTargetTemp()) {
    temperature = simParams->getTargetTemp();
  } else {
    sprintf(painCave.errMsg, "target temperature must be set\n");
    painCave.isFatal = 1;
    simError();  
  }
 
  std::map<std::string, SDShape> uniqueStuntDoubles;
  
  for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
      for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(ii)) {
          if (uniqueStuntDoubles.find(integrableObject->getType()) ==  uniqueStuntDoubles.end()) {

            integrableObject->setPos(V3Zero);
            integrableObject->setA(identMat);
            if (integrableObject->isRigidBody()) {
                RigidBody* rb = static_cast<RigidBody*>(integrableObject);
                rb->updateAtoms();
            }

            SDShape tmp;
            tmp.shape = ShapeBuilder::createShape(integrableObject);
            tmp.sd = integrableObject;    
            uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(integrableObject->getType(), tmp));

          }
        }
  }


  std::map<std::string, SDShape>::iterator si;
  for (si = uniqueStuntDoubles.begin(); si != uniqueStuntDoubles.end(); ++si) {
      HydrodynamicsModel* model;
      Shape* shape = si->second.shape;
      StuntDouble* sd = si->second.sd;;
      if (args_info.model_given) {  
        model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info);
      } else if (shape->hasAnalyticalSolution()) {
        model = new AnalyticalModel(sd, info);
      } else {
        model = new BeadModel(sd, info);
      }

        model->init();
        
        std::ofstream ofs;
        std::stringstream outputBeads;
        outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz";
        ofs.open(outputBeads.str().c_str());        
        model->writeBeads(ofs);
        ofs.close();

        //if beads option is turned on, skip the calculation
        if (!args_info.beads_flag) {
          model->calcHydroProps(shape, viscosity, temperature);
          std::ofstream outputDiff;
          outputDiff.open(outputFilename.c_str());
          model->writeHydroProps(outputDiff);
          outputDiff.close();
        }
        
        delete model;
  }


  //MemoryUtils::deletePointers(shapes);
  delete info;
   
}

void registerHydrodynamicsModels() {
    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel"));
    HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel"));

}
Revision:	1725
Committed:	Sat May 26 18:13:43 2012 UTC (12 years, 11 months ago) by gezelter
File size:	7194 byte(s)
Log Message:	Individual ForceField classes have been removed (they were essentially all duplicates anyway). ForceField has moved to brains, and since only one force field is in play at any time, the ForceFieldFactory and Register methods have been removed.
#	User	Rev	Content
1	tim	891	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9	gezelter	1390	* 1. Redistributions of source code must retain the above copyright
10	tim	891	* notice, this list of conditions and the following disclaimer.
11			*
12	gezelter	1390	* 2. Redistributions in binary form must reproduce the above copyright
13	tim	891	* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31	gezelter	1390	*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	tim	891	*/
42
43			#include <iostream>
44			#include <fstream>
45			#include <string>
46
47			#include "applications/hydrodynamics/HydroCmd.h"
48			#include "applications/hydrodynamics/HydrodynamicsModel.hpp"
49			#include "applications/hydrodynamics/HydrodynamicsModelCreator.hpp"
50			#include "applications/hydrodynamics/HydrodynamicsModelFactory.hpp"
51	tim	906	#include "applications/hydrodynamics/AnalyticalModel.hpp"
52	tim	891	#include "applications/hydrodynamics/BeadModel.hpp"
53			#include "applications/hydrodynamics/RoughShell.hpp"
54	tim	906	#include "applications/hydrodynamics/ShapeBuilder.hpp"
55	tim	891	#include "brains/Register.hpp"
56			#include "brains/SimCreator.hpp"
57			#include "brains/SimInfo.hpp"
58	tim	906	#include "utils/StringUtils.hpp"
59			#include "utils/simError.h"
60			#include "utils/MemoryUtils.hpp"
61	gezelter	1390	using namespace OpenMD;
62	tim	891
63	tim	906	struct SDShape{
64			StuntDouble* sd;
65			Shape* shape;
66			};
67	tim	891	void registerHydrodynamicsModels();
68	gezelter	972	void writeHydroProps(std::ostream& os);
69	tim	891
70			int main(int argc, char* argv[]){
71	tim	907	registerHydrodynamicsModels();
72	tim	891
73			gengetopt_args_info args_info;
74			std::string dumpFileName;
75			std::string mdFileName;
76			std::string prefix;
77
78			//parse the command line option
79			if (cmdline_parser (argc, argv, &args_info) != 0) {
80			exit(1) ;
81			}
82
83			//get the dumpfile name and meta-data file name
84			if (args_info.input_given){
85			dumpFileName = args_info.input_arg;
86			} else {
87			std::cerr << "Does not have input file name" << std::endl;
88			exit(1);
89			}
90
91			if (args_info.output_given){
92			prefix = args_info.output_arg;
93			} else {
94	gezelter	1027	prefix = getPrefix(dumpFileName);
95	tim	891	}
96	tim	892	std::string outputFilename = prefix + ".diff";
97	tim	906
98	tim	891	//parse md file and set up the system
99			SimCreator creator;
100	gezelter	1027	SimInfo* info = creator.createSim(dumpFileName, true);
101	tim	891
102			SimInfo::MoleculeIterator mi;
103			Molecule* mol;
104	tim	892	Molecule::IntegrableObjectIterator ii;
105			StuntDouble* integrableObject;
106			Mat3x3d identMat;
107			identMat(0,0) = 1.0;
108			identMat(1,1) = 1.0;
109			identMat(2,2) = 1.0;
110	tim	906
111			Globals* simParams = info->getSimParams();
112	tim	977	RealType temperature;
113			RealType viscosity;
114	tim	906
115			if (simParams->haveViscosity()) {
116			viscosity = simParams->getViscosity();
117			} else {
118	tim	907	sprintf(painCave.errMsg, "viscosity must be set\n");
119	tim	906	painCave.isFatal = 1;
120			simError();
121			}
122
123			if (simParams->haveTargetTemp()) {
124			temperature = simParams->getTargetTemp();
125			} else {
126	gezelter	945	sprintf(painCave.errMsg, "target temperature must be set\n");
127	tim	906	painCave.isFatal = 1;
128			simError();
129			}
130	tim	892
131	tim	906	std::map<std::string, SDShape> uniqueStuntDoubles;
132	tim	892
133	tim	891	for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
134	tim	892	for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
135			integrableObject = mol->nextIntegrableObject(ii)) {
136			if (uniqueStuntDoubles.find(integrableObject->getType()) == uniqueStuntDoubles.end()) {
137	tim	906
138	tim	892	integrableObject->setPos(V3Zero);
139			integrableObject->setA(identMat);
140			if (integrableObject->isRigidBody()) {
141			RigidBody* rb = static_cast<RigidBody*>(integrableObject);
142			rb->updateAtoms();
143			}
144	tim	921
145			SDShape tmp;
146			tmp.shape = ShapeBuilder::createShape(integrableObject);
147			tmp.sd = integrableObject;
148			uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(integrableObject->getType(), tmp));
149
150	tim	892	}
151			}
152	tim	891	}
153	tim	892
154	tim	906
155
156			std::map<std::string, SDShape>::iterator si;
157	tim	907	for (si = uniqueStuntDoubles.begin(); si != uniqueStuntDoubles.end(); ++si) {
158	tim	906	HydrodynamicsModel* model;
159			Shape* shape = si->second.shape;
160			StuntDouble* sd = si->second.sd;;
161			if (args_info.model_given) {
162			model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info);
163			} else if (shape->hasAnalyticalSolution()) {
164			model = new AnalyticalModel(sd, info);
165			} else {
166			model = new BeadModel(sd, info);
167			}
168	tim	921
169			model->init();
170
171			std::ofstream ofs;
172			std::stringstream outputBeads;
173			outputBeads << prefix << "_" << model->getStuntDoubleName() << ".xyz";
174			ofs.open(outputBeads.str().c_str());
175			model->writeBeads(ofs);
176			ofs.close();
177
178			//if beads option is turned on, skip the calculation
179			if (!args_info.beads_flag) {
180	gezelter	972	model->calcHydroProps(shape, viscosity, temperature);
181	xsun	1177	std::ofstream outputDiff;
182			outputDiff.open(outputFilename.c_str());
183	gezelter	972	model->writeHydroProps(outputDiff);
184	xsun	1177	outputDiff.close();
185	tim	921	}
186	gezelter	972
187	tim	906	delete model;
188	tim	891	}
189	tim	906
190
191			//MemoryUtils::deletePointers(shapes);
192	tim	891	delete info;
193
194			}
195
196			void registerHydrodynamicsModels() {
197			HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<RoughShell>("RoughShell"));
198			HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<BeadModel>("BeadModel"));
199	tim	906	HydrodynamicsModelFactory::getInstance()->registerHydrodynamicsModel(new HydrodynamicsModelBuilder<AnalyticalModel>("AnalyticalModel"));
200	tim	891
201			}