Changeset 9a2c8c for pcp/src/ions.c

Timestamp:

Apr 21, 2008, 2:19:25 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

4782e78

Parents:

8dd187

git-author:

Frederik Heber <heber@…> (04/18/08 16:16:38)

git-committer:

Frederik Heber <heber@…> (04/21/08 14:19:25)

Message:

CorrectForces(): greater fixes

File:

• pcp/src/ions.c (modified) (4 diffs)

pcp/src/ions.c

@@ -733 +733 @@
 
 /** Shifts center of gravity of ion forces IonType::FIon.
- * First sums up all forces of movable ions to a "force temperature",
- * then reduces each force by this temp, so that all in all
+ * First sums up all forces of movable ions for net center of gravity force,
+ * then reduces each force by this temp over Ions#Max_TotalIons, so that all in all
  * the net force is 0.
  * \param *P Problem at hand
- * \sa CorrectVelocity
- * \note why is FTemp not divided by number of ions: is probably correct, but this
- *       function was not really meaningful from the beginning.
+ * \sa CorrectVelocity()
  */
 void CorrectForces(struct Problem *P)
@@ -746 +744 @@
   double *FIon;
   double FTemp[NDIM];
-  int is, ia, d;
+  int is, ia, d, NumIons = 0;
   for (d=0; d<NDIM; d++)
-    FTemp[d] = 0;
+    FTemp[d] = 0.;
   for (is=0; is < I->Max_Types; is++) {         // calculate force temperature for each type ...
     for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {    // .. and each ion of this type ...
       FIon = &I->I[is].FIon[NDIM*ia];
-      if (I->I[is].IMT[ia] == MoveIon) {                // .. if it's movable
+      if (I->I[is].IMT[ia] == MoveIon) {                // .. if it's influenced
+        NumIons++;
                                 for (d=0; d<NDIM; d++) 
                                   FTemp[d] += FIon[d];
@@ -758 +757 @@
     }
   }
+  //if (P->Files.MeOutMes != 1) return;
+  //  fprintf(stderr, "TotalForce before: %e\t %e\t %e\n", FTemp[0], FTemp[1], FTemp[2]);
   for (is=0; is < I->Max_Types; is++) {         // and then reduce each by this value
     for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {
@@ -763 +764 @@
       if (I->I[is].IMT[ia] == MoveIon) {
                                 for (d=0; d<NDIM; d++) 
-                                 FIon[d] -= FTemp[d];                                           // (?) why not -= FTemp[d]/I->Max_TotalIons
-      }
-    }
-  }
+                                 FIon[d] -= FTemp[d]/NumIons;
+      }
+    }
+  }
+/*  for (d=0; d<NDIM; d++)
+    FTemp[d] = 0.;
+  for (is=0; is < I->Max_Types; is++) {         // calculate force temperature for each type ...
+    for (ia=0; ia < I->I[is].Max_IonsOfType; ia++) {    // .. and each ion of this type ...
+      FIon = &I->I[is].FIon[NDIM*ia];
+      if (I->I[is].IMT[ia] == MoveIon) {                // .. if it's influenced
+                                for (d=0; d<NDIM; d++) 
+                                  FTemp[d] += FIon[d];
+      }
+    }
+  }
+  if (P->Files.MeOutMes != 1) return;
+    fprintf(stderr, "TotalForce after: %e\t %e\t %e\n", FTemp[0], FTemp[1], FTemp[2]);
+    */
 }