vmaxwelln.cpp

pic 模拟程序！面向对象

/*
 ====================================================================
 VMAXWELLN.CPP
 0.99 (KC, 01-14-96) birth of a new class
 ====================================================================
*/

#include "vmaxwelln.h"
#include "ptclgrp.h"
#include "boundary.h"

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <txc_streams.h>
using namespace std;

#ifndef DBL_MIN
#define DBL_MIN         1E-200
#endif

//--------------------------------------------------------------------
void MaxwellianFlux::CreateZeroDist()
{
	FluxArray = NULL;
}
//--------------------------------------------------------------------
//	Set vt to the specified vector.  Default units are m/s.

void MaxwellianFlux::setThermal(const Vector3& v, UNIT type)
{
	Maxwellian::setThermal(v,type);
}
//--------------------------------------------------------------------
//	Set vc to the specified vector.  Default units are m/s.

void MaxwellianFlux::setCutoff(const Vector3& vl,const Vector3& vu, UNIT type)
{
	vcl = vl;
	vcu = vu;

	if ((vcl.magnitude()==0) && (vcu.magnitude()==0)){
		cutoffFlag = 0;
	}
	else{
		// brut force method to pass the neg sign through the sqrt
		Vector3 sign1,sign2;
		for (int e=1; e<=3;e++){
			if (vcl.e(e) == fabs(vcl.e(e)))
				sign1.set(e,1);
			else{
				vcl.set(e,fabs(vcl.e(e)));
				sign1.set(e,-1);
			}
			if (vcu.e(e) == fabs(vcu.e(e)))
				sign2.set(e,1);
			else{
				vcu.set(e,fabs(vcu.e(e)));
				sign2.set(e,-1);
			}
		}

		cutoffFlag = 1;
		//	If the units are in eV, we must convert to m/s internally:
		if(type==EV) {
				vcl *= 2*fabs(ELECTRON_CHARGE)/species->get_m();
				vcl = vcl.ComponentSqrt();
				vcl *= sign1;
				
				vcu *= 2*fabs(ELECTRON_CHARGE)/species->get_m();
				vcu = vcu.ComponentSqrt();
				vcu *= sign2;
		}
                else //SK for MKS the negative velocity limits disappeared!
                {
                       vcl *= sign1;
                       vcu *= sign2;
                }
		vcl = vcl.jvDivide(vt);
		vcu = vcu.jvDivide(vt);
		Maxwellian::ArraySetUp();
		if (v0.magnitude() ==0){
			if (fabs(vcu.magnitude())>3.3)
				cout << "Upper cutoff larger than most simulation will be able to resolve." << endl;
			expvu2 = Vector3(exp(sqr(vcu.e1())), exp(sqr(vcu.e2())), exp(sqr(vcu.e3())));
			expvl2 = Vector3(exp(sqr(vcl.e1())), exp(sqr(vcl.e2())), exp(sqr(vcl.e3())));
			Vsqr = vcl.ComponentSqr()+vcu.ComponentSqr();
		}
		else if ((v0.magnitude())<(50*vt.magnitude()))
			ArraySetUp();
	}
}

//--------------------------------------------------------------------
//	Set v0 to the specified vector.  Default units are m/s.

void MaxwellianFlux::setDrift(const Vector3& v, UNIT type)
{
    DriftFlag = 0; //default flag, used if length of V is 0
                   // David, please verify that it is correct m.g.
	Maxwellian::setDrift(v,type);
	if (v0.magnitude()==0)
		return;
	else{
		DriftFlag = 1;
		if ((v0.magnitude())<=(50*vt.magnitude()))
			ArraySetUp();
	}
}

//--------------------------------------------------------------------
//	Return a new velocity in a distribution with vt, vc, and v0 for 
// 	Maxwellian flux

Vector3 MaxwellianFlux::get_V(const Vector3& wall_normal)
{
	Vector3 v;
	Scalar vmag = 0;
	if ((v0.magnitude())>=(50*vt.magnitude()))
		return (Maxwellian::get_V());
	if ((cutoffFlag)&&(DriftFlag==0)){
		v = Maxwellian::get_V();
		Scalar f = frand();
		if ((wall_normal.e1()==1)||(wall_normal.e1()==-1))
			vmag = vt.e1()*sqrt(Vsqr.e1()+fabs(log(fabs(f*expvl2.e1()+(1-f)*expvu2.e1()))));
		else
			vmag = vt.e2()*sqrt(Vsqr.e2()+fabs(log(fabs(f*expvl2.e2()+(1-f)*expvu2.e2()))));
		if (wall_normal.e1()==1)
			v.set_e1(vmag);
		else if (wall_normal.e1()==-1)
			v.set_e1(-vmag);
		else if (wall_normal.e2()==1)
			v.set_e2(vmag);
		else if (wall_normal.e2()==-1)
			v.set_e2(-vmag);
	}
	else if ((DriftFlag)&&(cutoffFlag)){
		v = Maxwellian::get_V();
		Scalar temp1 = (nvel-1)*frand2();
		if ((wall_normal.e1()==1)||(wall_normal.e1()==-1))
			vmag =vt.e1()*((1-temp1+(int)temp1)*FluxArray[(int)temp1].e1()+(temp1-(int)temp1)*FluxArray[(int)temp1+1].e1());
		else
			vmag =vt.e2()*((1-temp1+(int)temp1)*FluxArray[(int)temp1].e2()+(temp1-(int)temp1)*FluxArray[(int)temp1+1].e2());
		if (wall_normal.e1()==1)
			v.set_e1(vmag);
		else if (wall_normal.e1()==-1)
			v.set_e1(-vmag);
		else if (wall_normal.e2()==1)
			v.set_e2(vmag);
		else if (wall_normal.e2()==-1)
			v.set_e2(-vmag);
	}
	else {
		if (DriftFlag){
			Scalar temp1 = (nvel-1)*frand();
			temp1 = (nvel-1)*frand();
			if ((wall_normal.e1()==1)||(wall_normal.e1()==-1))
				vmag =((1-temp1+(int)temp1)*FluxArray[(int)temp1].e1()+(temp1-(int)temp1)*FluxArray[(int)temp1+1].e1());
			else
				vmag =((1-temp1+(int)temp1)*FluxArray[(int)temp1].e2()+(temp1-(int)temp1)*FluxArray[(int)temp1+1].e2()); 
		}	
		else
			vmag= sqrt(2.0*fabs(log(fabs(frand()+1e-7))));
		
		if (wall_normal.e1()==1)
			v= (vt.jvMult(Vector3(vmag,normal(),normal())));
		else if (wall_normal.e1()==-1)
			v=(vt.jvMult(Vector3(-vmag,normal(),normal())));
		else if (wall_normal.e2()==1)
			v=(vt.jvMult(Vector3(normal(),vmag,normal())));	
		else if (wall_normal.e2()==-1)
			v=(vt.jvMult(Vector3(normal(),-vmag,normal())));
	}
	return v;
	
}

//--------------------------------------------------------------------
//	Return a new velocity in a distribution with vt, vc, and v0 for 
// 	Maxwellian flux

Vector3 MaxwellianFlux::get_U(const Vector3& wall_normal)
{
	Vector3 v = get_V(wall_normal);
	return (gamma0*v);
	
}

void MaxwellianFlux::ArraySetUp()
{
	nvel = 1000;
	FluxArray = new Vector3[nvel];
	Scalar v0local, vcllocal, vculocal;
	Scalar t1,t2,t3;
	Scalar dvi,dv,f,vtemp;
	Scalar flower,fupper,vlower,vupper,fmid;

	if ((vcl.magnitude() ==0)&&(vcu.magnitude()==0))
		{
			FluxArray[0] = Vector3(0,0,0);
			for (int e=1; e<=3; e++){
				v0local = v0.e(e)/vt.e(e);
				dvi = 3.3/((Scalar)(nvel-1));
				t2 = sqrt(M_PI)*v0local;
				t1 = exp(-sqr(v0local))+t2*erf(v0local);
				
				dv = dvi;
				fmid = 0;
				f = 0;
				vlower = 0;
				vtemp = vlower-v0local;
				flower = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t2*erf(vtemp);
				vupper = vlower + dv;
				vtemp = vupper-v0local;
				fupper = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t2*erf(vtemp);
				for (int n=1; n<nvel; n++){
					f = (Scalar)n/(Scalar)(nvel-1);
					vtemp = vlower-v0local;
					flower = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t2*erf(vtemp);
					vtemp = vupper-v0local;
					fupper = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t2*erf(vtemp);
					while ((flower*fupper)>0.0){
						dv*=2;
						vupper = vlower + dv;
						vtemp = vupper-v0local;
						fupper = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t2*erf(vtemp);
					}
					while (dv>1e-8){
						dv /=2;
						vtemp = vlower+dv-v0local;
						fmid = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t2*erf(vtemp);
						if (fmid<=0) 
							vlower+=dv;
						else
							vupper-=dv;
					}
					FluxArray[n].set(e,vlower); 
					dv = dvi;
				}
			}
		}
	else
		{
			for (int e=1; e<=3; e++){
				v0local = v0.e(e)/vt.e(e);
				vcllocal = vcl.e1()+v0local;
				vculocal = vcu.e1()+v0local;
				
				if (vcllocal<0)
					vcllocal=0;
				if (vculocal<0)
					vculocal=0;
				
				FluxArray[0].set(e,vcllocal);
				FluxArray[nvel-1].set(e,vculocal);
				dvi = (vculocal-vcllocal)/((Scalar)(nvel-1));
				t3 = sqrt(M_PI)*v0local;
				t1 = exp(-sqr(vcllocal-v0local))-t3*erf(vcllocal-v0local);
				t2 = -exp(-sqr(vculocal-v0local))+t3*erf(vculocal-v0local);
				
				dv = dvi;
				fmid = 0;
				vlower = vcllocal;
				f = 0;
				vtemp = vlower-v0local;
				flower = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t3*erf(vtemp);
				vupper = vlower + dv;
				vtemp = vupper-v0local;
				fupper = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t3*erf(vtemp);
				// for the 1 direction
				for (int n=1; n<nvel; n++){
					f = (Scalar)n/((Scalar)(nvel-1));
					vtemp = vlower-v0local;
					flower = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t3*erf(vtemp);
					vtemp = vupper-v0local;
					fupper = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t3*erf(vtemp);
					while ((flower*fupper)>0.0){
						dv*=2;
						vupper = vlower + dv;
						vtemp = vupper-v0local;
						fupper = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t3*erf(vtemp);
					}
					while (dv>1e-8){
						dv /=2;
						vtemp = vlower+dv-v0local;
						fmid = (1-f)*t1-f*t2-exp(-sqr(vtemp))+t3*erf(vtemp);
						if (fmid<=0) 
							vlower+=dv;
						else
							vupper-=dv;
					}
					FluxArray[n].set(e,vlower);
					dv=dvi;
				}
			}
		}		
}