Discussion Forum

In case it's interesting to others, the resolution to this problem is that COMSOL does not pass velocities for all particles to the external function, it only passes velocities for those particles which had a collision during the last timestep. If you pass cpt.idx along with your velocities, COMSOL is smart enough to only pass the indices that are relevant and all the vectors it passes are ordered such that the subset of cpt.idx properly indexes the other vectors.

The following C++ code will take some arguments from the charged particle tracing module and randomize one velocity component. It writes all the input and output variables to a file called variableStates so you can see what's going on under the hood.

Edit: sorry about the formatting - I can't figure out how to get the forum to keep spaces or tabs at the beginning of lines.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifdef _MSC_VER
#define EXPORT __declspec(dllexport)
#else
#define EXPORT
#endif

#ifdef __cplusplus
extern "C"
{
#endif

static const char *error = NULL;

EXPORT int init(const char *str) {
/*open a file and write a line indicating that the function init
has been called */
FILE *ofp;
char outputFilename[] = "variableStates";
ofp = fopen(outputFilename, "a");
fprintf(ofp, "%s\n", "--- init Called ---");

//done with that, close the file
fclose(ofp);
return 1;

}

EXPORT const char * getLastError() {
/*open a file and write a line indicating that the function
getLastError has been called */
FILE *ofp;
char outputFilename[] = "variableStates";
ofp = fopen(outputFilename, "a");
fprintf(ofp, "%s\n", "--- getLastError Called ---");

//done with that, close the file
fclose(ofp);
return error;

}

double sq(const double x){
//this is a function that squares its argument
return x*x;
}

EXPORT int eval(const char *func, int nArgs, const double **inReal, const double **inImag, int blockSize, double *outReal, double *outImag) {
/* required arguments:
cartesian velocity (cpt.vx, cpt.vy, or cpt.vz)
cpt.pidx
cpt.Nc_ecX where X is a sequential index for each elastic
collision force node in the Charged Particle Tracing module */

int i, idx, thisIdx, r;
double x;

//check if COMSOL asked for the right function
if (strcmp("ext12", func) == 0) {

/*check if COMSOL passed the right number of arguments
and return if it didn't */
if (nArgs != 3) {
error = "Two arguments expected";
return 0;
}

/*open a file and write a line indicating that the function
eval has been called, then dump some input variables */
FILE *ofp;
char outputFilename[] = "variableStates";
ofp = fopen(outputFilename, "a");
fprintf(ofp, "%s\n", "--- eval Called ---");
fprintf(ofp, "nArgs\t\t = %d\nblockSize\t = %d", nArgs, blockSize);

// write the particle indices to file
fprintf(ofp, "\nindices\t\t = ");
for (i = 0; i < blockSize; i++) {
fprintf(ofp, "%.0f\t\t", inReal[1][i]);
}

// write the numbers of collisions to file
fprintf(ofp, "\ncollisions\t = ");
for (i = 0; i < blockSize; i++) {
fprintf(ofp, "%.0f\t\t", inReal[2][i]);
}

// write the input velocities to file
fprintf(ofp, "\ninReal[0]\t = ");
for (i = 0; i < blockSize; i++) {
fprintf(ofp, "%f ", inReal[0][i]);
}

//randomize the velocity of each colliding particle
for (i = 0; i < blockSize; i++) {

/*set the current particle index and the number of
collisions the current particle has seen */
idx = inReal[1][i];

/*seed the random number generator with a function of
the particle index and the number of collisions
so that each particle gets a unique sequence of random
numbers each time it collides */
thisIdx = idx + inReal[2][i];
srand((thisIdx+8971)*(thisIdx+8971));

/*set the output velocity to a random number between
0 and 10 */
r = rand()%10000;
outReal[i] = r/1000.0;
}

// write the output velocities to file
fprintf(ofp, "\noutReal\t\t = ");
for (i = 0; i < blockSize; i++) {
fprintf(ofp, "%f ", outReal[i]);
}

//add a blank line at the end of the file and close it
fprintf(ofp, "\n");
fclose(ofp);

return 1;
}

//if COMSOL called the wrong function, say so and return
else {
error = "Unknown function";
return 0;
}

}

#ifdef __cplusplus
} // __cplusplus defined.
#endif