<commit_before>/*===- Runge_Kutta.cpp - libSimulation -========================================
*
*                                  DEMON
*
* This file is distributed under the BSD Open Source License. See LICENSE.TXT 
* for details.
*
*===-----------------------------------------------------------------------===*/

#include "Runge_Kutta.h"
#include <cmath>
#include <limits>
#include "VectorCompatibility.h"
#include "PositionVelocityCacheOperator.h"
#include "ChargeOperator.h"
#include "FieldPotentialOperator.h"

using namespace std;

Runge_Kutta::Runge_Kutta(Cloud * const myCloud, Force **forces, const double timeStep, const force_index forcesSize, const double startTime)
: cloud(myCloud), theForce(forces), numForces(forcesSize), init_dt(timeStep), currentTime(startTime), 
numOperators(3), operations(new Operator*[numOperators])
{
	// Operators are order dependent.
	operations[0] = new PositionVelocityCacheOperator(cloud);
	operations[1] = new FieldPotentialOperator(cloud);
	operations[2] = new ChargeOperator(cloud);
}

Runge_Kutta::~Runge_Kutta()
{
	for (operator_index i = 0; i < numOperators; i++)
		delete operations[i];
	delete[] operations;
}

// 4th order Runge-Kutta algorithm:
void Runge_Kutta::moveParticles(const double endTime)
{
	// create vector constants:
	const __m128d v2 = _mm_set1_pd(2.0);
	const __m128d v6 = _mm_set1_pd(6.0);
    
	while (currentTime < endTime)
	{
		const double dt = modifyTimeStep(0, 2, 1.0e-4, init_dt); // implement dynamic timstep (if necessary):
			// Second argument must be 2 more than the first.
		const __m128d vdt = _mm_set1_pd(dt); // store timestep as vector const

		const double const1 = cloud->chargeConst1;
		const double const2 = cloud->chargeConst2;
		const double particleRadius = cloud->particleRadius;
		const __m128d qConst1 = _mm_set_pd(const1, const1);
		const __m128d qConst2 = _mm_set_pd(const2, const2);
		const __m128d partRadius = _mm_set_pd(particleRadius, particleRadius);
        
		operate1(currentTime);
		force1(currentTime); // compute net force1
		for (cloud_index i = 0, numParticles = cloud->n; i < numParticles; i += 2) // calculate k1 and l1 for entire cloud
		{
			const __m128d vmass = _mm_load_pd(cloud->mass + i); // load ith and (i+1)th mass into vector

			// assign pointers:
			double * const pFx = cloud->forceX + i;
			double * const pFy = cloud->forceY + i;
			double * const pPhi = cloud->phi + i;
			double * const pEx = cloud->Ex + i;
			double * const pEy = cloud->Ey + i;
           
			// calculate ith and (i+1)th tidbits: 
			_mm_store_pd(cloud->k1 + i, vdt*_mm_load_pd(pFx)/vmass); // velocityX tidbit
			_mm_store_pd(cloud->l1 + i, vdt*cloud->getVx1_pd(i)); // positionX tidbit
			_mm_store_pd(cloud->m1 + i, vdt*_mm_load_pd(pFy)/vmass); // velocityY tidbit
			_mm_store_pd(cloud->n1 + i, vdt*cloud->getVy1_pd(i)); // positionY tidbit
			//_mm_store_pd(cloud->q1 + i, -vdt*(qConst1*_mm_load_pd(cloud->charge + i) +
			//	qConst2*partRadius*_mm_load_pd(cloud->phi + i)));
			_mm_store_pd(cloud->q1 + i, _mm_load_pd(cloud->charge + i));
				// see "Instability of Dust Acoustic Waves in an Accelerating Dusy Plasma"
				// <NAME>, <NAME>, <NAME>
				// Physica Scripta, Vol. 67, 354-356, 2003

			// reset to zero:
			_mm_store_pd(pFx, _mm_setzero_pd());
			_mm_store_pd(pFy, _mm_setzero_pd());
			_mm_store_pd(pPhi, _mm_setzero_pd());
			_mm_store_pd(pEx, _mm_setzero_pd());
			_mm_store_pd(pEy, _mm_setzero_pd());
			
		}
        
		operate2(currentTime + dt/2.0);
		force2(currentTime + dt/2.0); // compute net force2
		for (cloud_index i = 0, numParticles = cloud->n; i < numParticles; i += 2) // calculate k2 and l2 for entire cloud
		{
			const __m128d vmass = _mm_load_pd(cloud->mass + i); // load ith and (i+1)th mass

			// assign pointers:
			double * const pFx = cloud->forceX + i;
			double * const pFy = cloud->forceY + i;
			double * const pPhi = cloud->phi + i;
			double * const pEx = cloud->Ex + i;
			double * const pEy = cloud->Ey + i;

			// calculate ith and (i+1)th tidbits: 
			_mm_store_pd(cloud->k2 + i, vdt*_mm_load_pd(pFx)/vmass); // velocityX tidbit
			_mm_store_pd(cloud->l2 + i, vdt*cloud->getVx2_pd(i)); // positionX tidbit
			_mm_store_pd(cloud->m2 + i, vdt*_mm_load_pd(pFy)/vmass); // velocityY tidbit
			_mm_store_pd(cloud->n2 + i, vdt*cloud->getVy2_pd(i)); // positionY tidbit
			//_mm_store_pd(cloud->q2 + i, -vdt*(qConst1*cloud->getq1_pd(i) +
			//	qConst2*partRadius*_mm_load_pd(cloud->phi + i)));
			_mm_store_pd(cloud->q2 + i, _mm_load_pd(cloud->charge + i));

			// reset to zero:
			_mm_store_pd(pFx, _mm_setzero_pd());
			_mm_store_pd(pFy, _mm_setzero_pd());
			_mm_store_pd(pPhi, _mm_setzero_pd());
			_mm_store_pd(pEx, _mm_setzero_pd());
			_mm_store_pd(pEy, _mm_setzero_pd());
		}

		operate3(currentTime + dt/2.0);
		force3(currentTime + dt/2.0); // compute net force3
		for (cloud_index i = 0, numParticles = cloud->n; i < numParticles; i += 2) // calculate k3 and l3 for entire cloud
		{
			const __m128d vmass = _mm_load_pd(cloud->mass + i); // load ith and (i+1)th mass

			// assign force pointers:
			double * const pFx = cloud->forceX + i;
			double * const pFy = cloud->forceY + i;
			double * const pPhi = cloud->phi + i;
			double * const pEx = cloud->Ex + i;
			double * const pEy = cloud->Ey + i;

			// calculate ith and (i+1)th tibits: 
			_mm_store_pd(cloud->k3 + i, vdt*_mm_load_pd(pFx)/vmass); // velocityX tidbit
			_mm_store_pd(cloud->l3 + i, vdt*cloud->getVx3_pd(i)); // positionX tidbit
			_mm_store_pd(cloud->m3 + i, vdt*_mm_load_pd(pFy)/vmass); // velocityY tidbit
			_mm_store_pd(cloud->n3 + i, vdt*cloud->getVy3_pd(i)); // positionY tidbit
			//_mm_store_pd(cloud->q3 + i, -vdt*(qConst1*cloud->getq2_pd(i) +
			//	qConst2*partRadius*_mm_load_pd(cloud->phi + i)));
			_mm_store_pd(cloud->q3 + i, _mm_load_pd(cloud->charge + i));

			// reset forces to zero:
			_mm_store_pd(pFx, _mm_setzero_pd());
			_mm_store_pd(pFy, _mm_setzero_pd());
			_mm_store_pd(pPhi, _mm_setzero_pd());
			_mm_store_pd(pEx, _mm_setzero_pd());
			_mm_store_pd(pEy, _mm_setzero_pd());
		}
        
		operate4(currentTime + dt/2.0);
		force4(currentTime + dt); // compute net force4
		for (cloud_index i = 0, numParticles = cloud->n; i < numParticles; i += 2) // calculate k4 and l4 for entire cloud
		{
			const __m128d vmass = _mm_load_pd(cloud->mass + i); // load ith and (i+1)th mass

			// assign force pointers:
			double * const pFx = cloud->forceX + i;
			double * const pFy = cloud->forceY + i;
			double * const pPhi = cloud->phi + i;
			double * const pEx = cloud->Ex + i;
			double * const pEy = cloud->Ey + i;
            
			_mm_store_pd(cloud->k4 + i, vdt*_mm_load_pd(pFx)/vmass); // velocityX tidbit
			_mm_store_pd(cloud->l4 + i, vdt*cloud->getVx4_pd(i)); // positionX tidbit
			_mm_store_pd(cloud->m4 + i, vdt*_mm_load_pd(pFy)/vmass); // velocityY tidbit
			_mm_store_pd(cloud->n4 + i, vdt*cloud->getVy4_pd(i)); // positionY tidbit
			//_mm_store_pd(cloud->q4 + i, -vdt*(qConst1*cloud->getq3_pd(i) +
			//	qConst2*partRadius*_mm_load_pd(cloud->phi + i)));
			_mm_store_pd(cloud->q4 + i, _mm_load_pd(cloud->charge + i));

			// reset forces to zero:
			_mm_store_pd(pFx, _mm_setzero_pd());
			_mm_store_pd(pFy, _mm_setzero_pd());
			_mm_store_pd(pPhi, _mm_setzero_pd());
			_mm_store_pd(pEx, _mm_setzero_pd());
			_mm_store_pd(pEy, _mm_setzero_pd());
		}<commit_msg>Place dynamic charge calculation is a preprocessor define. This allows the charge calculation to be turned on and off by placing -D CHARGE in the CXXFLAGS of the makefile.include<commit_after>/*===- Runge_Kutta.cpp - libSimulation -========================================
*
*                                  DEMON
*
* This file is distributed under the BSD Open Source License. See LICENSE.TXT 
* for details.
*
*===-----------------------------------------------------------------------===*/

#include "Runge_Kutta.h"
#include <cmath>
#include <limits>
#include "VectorCompatibility.h"
#include "PositionVelocityCacheOperator.h"
#include "ChargeOperator.h"
#include "FieldPotentialOperator.h"

using namespace std;

Runge_Kutta::Runge_Kutta(Cloud * const myCloud, Force **forces, const double timeStep, const force_index forcesSize, const double startTime)
: cloud(myCloud), theForce(forces), numForces(forcesSize), init_dt(timeStep), currentTime(startTime), 
numOperators(3), operations(new Operator*[numOperators])
{
	// Operators are order dependent.
	operations[0] = new PositionVelocityCacheOperator(cloud);
	operations[1] = new FieldPotentialOperator(cloud);
	operations[2] = new ChargeOperator(cloud);
}

Runge_Kutta::~Runge_Kutta()
{
	for (operator_index i = 0; i < numOperators; i++)
		delete operations[i];
	delete[] operations;
}

// 4th order Runge-Kutta algorithm:
void Runge_Kutta::moveParticles(const double endTime)
{
	// create vector constants:
	const __m128d v2 = _mm_set1_pd(2.0);
	const __m128d v6 = _mm_set1_pd(6.0);
    
	while (currentTime < endTime)
	{
		const double dt = modifyTimeStep(0, 2, 1.0e-4, init_dt); // implement dynamic timstep (if necessary):
			// Second argument must be 2 more than the first.
		const __m128d vdt = _mm_set1_pd(dt); // store timestep as vector const
#ifdef CHARGE
		const __m128d qConst1 = _mm_set1_pd(cloud->chargeConst1);
		const __m128d qConst2 = _mm_set1_pd(cloud->chargeConst2);
		const __m128d partRadius = _mm_set1_pd(cloud->particleRadius);
#endif
		
		operate1(currentTime);
		force1(currentTime); // compute net force1
		for (cloud_index i = 0, numParticles = cloud->n; i < numParticles; i += 2) // calculate k1 and l1 for entire cloud
		{
			const __m128d vmass = _mm_load_pd(cloud->mass + i); // load ith and (i+1)th mass into vector

			// assign pointers:
			double * const pFx = cloud->forceX + i;
			double * const pFy = cloud->forceY + i;
			double * const pPhi = cloud->phi + i;
			double * const pEx = cloud->Ex + i;
			double * const pEy = cloud->Ey + i;
           
			// calculate ith and (i+1)th tidbits: 
			_mm_store_pd(cloud->k1 + i, vdt*_mm_load_pd(pFx)/vmass); // velocityX tidbit
			_mm_store_pd(cloud->l1 + i, vdt*cloud->getVx1_pd(i)); // positionX tidbit
			_mm_store_pd(cloud->m1 + i, vdt*_mm_load_pd(pFy)/vmass); // velocityY tidbit
			_mm_store_pd(cloud->n1 + i, vdt*cloud->getVy1_pd(i)); // positionY tidbit
#ifdef CHARGE
			_mm_store_pd(cloud->q1 + i, -vdt*(qConst1*_mm_load_pd(cloud->charge + i) +
				qConst2*partRadius*_mm_load_pd(cloud->phi + i)));
#else
			_mm_store_pd(cloud->q1 + i, _mm_load_pd(cloud->charge + i));
#endif
				// see "Instability of Dust Acoustic Waves in an Accelerating Dusy Plasma"
				// <NAME>, <NAME>, <NAME>
				// Physica Scripta, Vol. 67, 354-356, 2003

			// reset to zero:
			_mm_store_pd(pFx, _mm_setzero_pd());
			_mm_store_pd(pFy, _mm_setzero_pd());
			_mm_store_pd(pPhi, _mm_setzero_pd());
			_mm_store_pd(pEx, _mm_setzero_pd());
			_mm_store_pd(pEy, _mm_setzero_pd());
			
		}
        
		operate2(currentTime + dt/2.0);
		force2(currentTime + dt/2.0); // compute net force2
		for (cloud_index i = 0, numParticles = cloud->n; i < numParticles; i += 2) // calculate k2 and l2 for entire cloud
		{
			const __m128d vmass = _mm_load_pd(cloud->mass + i); // load ith and (i+1)th mass

			// assign pointers:
			double * const pFx = cloud->forceX + i;
			double * const pFy = cloud->forceY + i;
			double * const pPhi = cloud->phi + i;
			double * const pEx = cloud->Ex + i;
			double * const pEy = cloud->Ey + i;

			// calculate ith and (i+1)th tidbits: 
			_mm_store_pd(cloud->k2 + i, vdt*_mm_load_pd(pFx)/vmass); // velocityX tidbit
			_mm_store_pd(cloud->l2 + i, vdt*cloud->getVx2_pd(i)); // positionX tidbit
			_mm_store_pd(cloud->m2 + i, vdt*_mm_load_pd(pFy)/vmass); // velocityY tidbit
			_mm_store_pd(cloud->n2 + i, vdt*cloud->getVy2_pd(i)); // positionY tidbit
			//_mm_store_pd(cloud->q2 + i, -vdt*(qConst1*cloud->getq1_pd(i) +
			//	qConst2*partRadius*_mm_load_pd(cloud->phi + i)));
			_mm_store_pd(cloud->q2 + i, _mm_load_pd(cloud->charge + i));

			// reset to zero:
			_mm_store_pd(pFx, _mm_setzero_pd());
			_mm_store_pd(pFy, _mm_setzero_pd());
			_mm_store_pd(pPhi, _mm_setzero_pd());
			_mm_store_pd(pEx, _mm_setzero_pd());
			_mm_store_pd(pEy, _mm_setzero_pd());
		}

		operate3(currentTime + dt/2.0);
		force3(currentTime + dt/2.0); // compute net force3
		for (cloud_index i = 0, numParticles = cloud->n; i < numParticles; i += 2) // calculate k3 and l3 for entire cloud
		{
			const __m128d vmass = _mm_load_pd(cloud->mass + i); // load ith and (i+1)th mass

			// assign force pointers:
			double * const pFx = cloud->forceX + i;
			double * const pFy = cloud->forceY + i;
			double * const pPhi = cloud->phi + i;
			double * const pEx = cloud->Ex + i;
			double * const pEy = cloud->Ey + i;

			// calculate ith and (i+1)th tibits: 
			_mm_store_pd(cloud->k3 + i, vdt*_mm_load_pd(pFx)/vmass); // velocityX tidbit
			_mm_store_pd(cloud->l3 + i, vdt*cloud->getVx3_pd(i)); // positionX tidbit
			_mm_store_pd(cloud->m3 + i, vdt*_mm_load_pd(pFy)/vmass); // velocityY tidbit
			_mm_store_pd(cloud->n3 + i, vdt*cloud->getVy3_pd(i)); // positionY tidbit
#ifdef CHARGE
			_mm_store_pd(cloud->q3 + i, -vdt*(qConst1*cloud->getq2_pd(i) +
				qConst2*partRadius*_mm_load_pd(cloud->phi + i)));
#else
			_mm_store_pd(cloud->q3 + i, _mm_load_pd(cloud->charge + i));
#endif
			
			// reset forces to zero:
			_mm_store_pd(pFx, _mm_setzero_pd());
			_mm_store_pd(pFy, _mm_setzero_pd());
			_mm_store_pd(pPhi, _mm_setzero_pd());
			_mm_store_pd(pEx, _mm_setzero_pd());
			_mm_store_pd(pEy, _mm_setzero_pd());
		}
        
		operate4(currentTime + dt/2.0);
		force4(currentTime + dt); // compute net force4
		for (cloud_index i = 0, numParticles = cloud->n; i < numParticles; i += 2) // calculate k4 and l4 for entire cloud
		{
			const __m128d vmass = _mm_load_pd(cloud->mass + i); // load ith and (i+1)th mass

			// assign force pointers:
			double * const pFx = cloud->forceX + i;
			double * const pFy = cloud->forceY + i;
			double * const pPhi = cloud->phi + i;
			double * const pEx = cloud->Ex + i;
			double * const pEy = cloud->Ey + i;
            
			_mm_store_pd(cloud->k4 + i, vdt*_mm_load_pd(pFx)/vmass); // velocityX tidbit
			_mm_store_pd(cloud->l4 + i, vdt*cloud->getVx4_pd(i)); // positionX tidbit
			_mm_store_pd(cloud->m4 + i, vdt*_mm_load_pd(pFy)/vmass); // velocityY tidbit
			_mm_store_pd(cloud->n4 + i, vdt*cloud->getVy4_pd(i)); // positionY tidbit
#ifdef CHARGE
			_mm_store_pd(cloud->q4 + i, -vdt*(qConst1*cloud->getq3_pd(i) +
				qConst2*partRadius*_mm_load_pd(cloud->phi + i)));
#else
			_mm_store_pd(cloud->q4 + i, _mm_load_pd(cloud->charge + i));
#endif
			
			// reset forces to zero:
			_mm_store_pd(pFx, _mm_setzero_pd());
			_mm_store_pd(pFy, _mm_setzero_pd());
			_mm_store_pd(pPhi, _mm_setzero_pd());
			_mm_store_pd(pEx, _mm_setzero_pd());
			_mm_store_pd(pEy, _mm_setzero_pd());
		}
{'fraction_non_alphanumeric': 0.1566786027971083, 'fraction_numerical': 0.02155259779744612, 'mean_word_length': 3.65764631843927, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 163, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1404263', 'n_tokens_mistral': 6412, 'n_tokens_neox': 5820, 'n_words': 1254}