C++ GridInfo::R方法代码示例

//################ Droplet on attractive wall ####################
int main(int argc, char** argv)
{	initSystem(argc, argv);

	GridInfo gInfo;
	gInfo.S = vector3<int>(64, 64, 64); double hGrid=1.0;
	//gInfo.S = vector3<int>(128, 128, 128); double hGrid=0.5;
	gInfo.R = Diag(gInfo.S * hGrid);
	gInfo.initialize();

	double T = 298*Kelvin;
	FluidComponent component(FluidComponent::H2O, T, FluidComponent::ScalarEOS);
	component.s2quadType = QuadOctahedron;
	component.Nnorm = 270;
	
	PreconditionedFluidMixture fluidMixture(gInfo, T, 1.0);
	component.addToFluidMixture(&fluidMixture);
	double p = 1.01325*Bar;
	logPrintf("pV = %le\n", p*gInfo.detR);
	fluidMixture.initialize(p);

	#define geomName "AttractiveWall-3bohr-3.0kT/drop_plane"

	bool loadState=false;
	const char stateFilename[] = "TestFixedN/" geomName "_state.bin";

	//Initialize potential: planar wall with attractive well:
	nullToZero(component.idealGas->V, gInfo);
	double xWall = 0.5*gInfo.R(0,0)-21.0;
	double dxWall = 2.0;
	applyFunc_r(gInfo, initAttractiveWall, xWall, dxWall, 100*T, 3.*T, component.idealGas->V[0]->data());

	if(loadState)
		fluidMixture.loadState(stateFilename);
	else
	{	//Initialize state biased towards center of cell
		const double muSet=-5.0;
		const double Rdroplet = 22.0; //guess droplet size:
		nullToZero(fluidMixture.state, gInfo, fluidMixture.get_nIndep());
		applyFunc_r(gInfo, initSphere, gInfo.R*vector3<>(0.5,0.5,0.5), Rdroplet, 0.0, muSet, fluidMixture.state[0]->data());

		RealKernel gauss(gInfo); initGaussianKernel(gauss, 2.0);
		fluidMixture.state[0] = I(gauss*J(fluidMixture.state[0]));

		ScalarField wallMask(ScalarFieldData::alloc(gInfo));
		applyFunc_r(gInfo, initAttractiveWall, xWall, 2*dxWall, 0.0, 1.0, wallMask->data());
		fluidMixture.state[0] *= (wallMask+1.0);
	}

	MinimizeParams mp;
	mp.fpLog = globalLog;
	mp.nDim = 2*gInfo.nr;
	mp.nIterations=10;
	mp.knormThreshold=1e-11;
	mp.dirUpdateScheme = MinimizeParams::FletcherReeves;
	//mp.dirUpdateScheme = MinimizeParams::SteepestDescent;
	//mp.updateTestStepSize = false;
	mp.fdTest = !loadState;

	int sysRet=system("mkdir -p TestFixedN/" geomName "_img");
	if(sysRet) { logPrintf("Error making image directory\n"); mpiUtil->exit(sysRet); }

	for(int loopCount=0; loopCount<100; loopCount++)
	{
		ScalarFieldArray N;
		TIME("getOmega calculation (with gradient)", globalLog,
			double omega = fluidMixture.getFreeEnergy(FluidMixture::Outputs(&N));
			if(std::isnan(omega)) break; //Don't save state after it has become nan
		);
		logPrintf("Ntot = %lf\n", gInfo.dV*sum(N[0]));

		logPrintf("Saving state:\n");
		fluidMixture.saveState(stateFilename);
		saveDX(N[0], "TestFixedN/" geomName "_nO");
		saveDX(N[1], "TestFixedN/" geomName "_nH");
		saveSphericalized(&N[0], 2, "TestFixedN/" geomName "_n.spherical", 0.25);
		//Invoke octave to create an image:
		FILE* pp = popen("octave -q", "w");
		fprintf(pp, "imwrite( waterSlice(\"TestFixedN/" geomName "_n%%s.bin\", [%d %d %d], 1, %d, 1e-2),", gInfo.S[0], gInfo.S[1], gInfo.S[2], gInfo.S[2]/2);
		fprintf(pp, " \"TestFixedN/" geomName "_img/img%04d.png\"); exit;\n", loopCount);
		fflush(pp); pclose(pp);

		logPrintf("Starting CG:\n");
		TIME("minimize", globalLog,
			fluidMixture.minimize(mp);
		);

int main(int argc, char** argv)
{	initSystem(argc, argv);

	//Parse command-line
	S2quadType quadType = QuadEuler;
	int nBeta = 12;
	if(argc > 1)
	{	if(!S2quadTypeMap.getEnum(argv[1], quadType))
			die("<quad> must be one of %s\n", S2quadTypeMap.optionList().c_str());
		if(quadType==QuadEuler)
		{	if(argc < 3) die("<nBeta> must be specified for Euler quadratures.\n")
			nBeta = atoi(argv[2]);
			if(nBeta <= 0) die("<nBeta> must be non-negative.")
		}
	}
	
	//Setup simulation grid:
	GridInfo gInfo;
	gInfo.S = vector3<int>(1, 1, 4096);
	const double hGrid = 0.0625;
	gInfo.R = Diag(hGrid * gInfo.S);
	gInfo.initialize();

	double T = 298*Kelvin;
	FluidComponent component(FluidComponent::H2O, T, FluidComponent::ScalarEOS);
	component.s2quadType = quadType;
	component.quad_nBeta = nBeta;
	component.representation = FluidComponent::Pomega;
	FluidMixture fluidMixture(gInfo, T);
	component.addToFluidMixture(&fluidMixture);
	double p = 1.01325*Bar;
	logPrintf("pV = %le\n", p*gInfo.detR);
	fluidMixture.initialize(p);

	//Initialize external potential (repel O from a cube)
	double Dfield = 1.0 * eV/Angstrom;
	const double zWall = 8.0 - 1e-3;
	const double& gridLength = gInfo.R(2,2);
	ScalarField phiApplied(ScalarFieldData::alloc(gInfo)), phiWall(ScalarFieldData::alloc(gInfo));
	applyFunc_r(gInfo, setPhi, phiApplied->data(), phiWall->data(), gridLength, Dfield, zWall);
	const double ZO = component.molecule.sites[0]->chargeKernel(0);
	component.idealGas->V[0] = ZO * phiApplied + phiWall;
	component.idealGas->V[1] = -0.5*ZO * phiApplied + phiWall;

	//----- Initialize state -----
	fluidMixture.initState(0.01);

	//----- FDtest and CG -----
	MinimizeParams mp;
	mp.fpLog = globalLog;
	mp.nDim = gInfo.nr * fluidMixture.get_nIndep();
	mp.energyLabel = "Phi";
	mp.nIterations=1500;
	mp.energyDiffThreshold=1e-16;
	
	fluidMixture.minimize(mp);
	
	//------ Outputs ---------
	ostringstream quadName;
	quadName << S2quadTypeMap.getString(quadType);
	if(quadType == QuadEuler) quadName << nBeta;
	
	ScalarFieldArray N;
	fluidMixture.getFreeEnergy(FluidMixture::Outputs(&N));
	FILE* fp = fopen((quadName.str()+".Nplanar").c_str(), "w");
	double* NOdata = N[0]->data();
	double* NHdata = N[1]->data();
	double nlInv = 1./component.idealGas->get_Nbulk();
	for(int i=0; i<gInfo.S[2]/2; i++)
		fprintf(fp, "%le\t%le\t%le\n", i*hGrid, nlInv*NOdata[i], 0.5*nlInv*NHdata[i]);
	fclose(fp);
	
	finalizeSystem();
	return 0;
}