Changeset 1681 for branches/newmole/source

Timestamp:

12/16/07 06:22:43 (11 months ago)

Author:

rjrw

Message:

Fixes to get dynamics running with conservation constraints applied.

ion_solver now renormalizes for /all/ cases, the dynamic sources to
the ladder are applied in the molecular network.

Small recodings to fold trimmed advective species into the highest
state available and reorder a logical test so the most likely failure
comes first (small short-circuit optimization).

Location:

branches/newmole/source

Files:

• dynamics.cpp (modified) (4 diffs)
• ion_solver.cpp (modified) (2 diffs)
• mole_newton_step.cpp (modified) (4 diffs)

branches/newmole/source/dynamics.cpp

@@ -774 +774 @@
                         {
                                 dynamics.Source[nelem][ion] = 0.;
+                                dynamics.Source[nelem][dense.IonHigh[nelem]] += 
+                                        UpstreamIon[nelem][ion] * scalingDensity() / timestep;
                         }
                 }
@@ -948 +950 @@
                                  Old_molecules[ipUpstream][mol]/Old_density[ipUpstream])*
                                 frac_next;
-                        /* Atoms and positive ions are already counted in upstreamElem */
+                        /* Atoms and positive ions are already counted in UpstreamElem */
                         if(COmole[mol]->n_nuclei != 1 || COmole[mol]->nElec < 0) 
                         {
@@ -1506 +1508 @@
                                         Oldi_ion = Old_xIonDense[ilast][nelem][ion];
                                 }
-                                dynamics.convergence_error += POW2(Oldi_ion/Oldi_density-struc.xIonDense[nelem][ion][i]/scalingZoneDensity(i)) /* *struc.dr[i] */;  
+                                dynamics.convergence_error += POW2((double)Oldi_ion/Oldi_density-struc.xIonDense[nelem][ion][i]/scalingZoneDensity(i)) /* *struc.dr[i] */;  
 
                                 /* >>chng 02 nov 11, add first error scale estimate from Robin */
-                                dynamics.error_scale1 += POW2(struc.xIonDense[nelem][ion][i]/scalingZoneDensity(i));                    
+                                //fprintf(ioQQQ,"%g %g %g\n",dynamics.error_scale1,
+                                //                              struc.xIonDense[nelem][ion][i],scalingZoneDensity(i));
+                                dynamics.error_scale1 += POW2((double)struc.xIonDense[nelem][ion][i]/scalingZoneDensity(i));                    
                         }
                 }
@@ -1524 +1528 @@
                                 Oldi_mol = Old_molecules[ilast][mol];
                         }
-                        dynamics.convergence_error += POW2(Oldi_mol/Oldi_density-struc.molecules[mol][i]/scalingZoneDensity(i)) /* *struc.dr[i] */;  
+                        dynamics.convergence_error += POW2((double)Oldi_mol/Oldi_density-struc.molecules[mol][i]/scalingZoneDensity(i)) /* *struc.dr[i] */;  
 
                         /* >>chng 02 nov 11, add first error scale estimate from Robin
                          * used to normalize the above convergence_error */
-                        dynamics.error_scale1 += POW2(struc.molecules[mol][i]/scalingZoneDensity(i));                   
+                        dynamics.error_scale1 += POW2((double)struc.molecules[mol][i]/scalingZoneDensity(i));                   
                 }
         }

branches/newmole/source/ion_solver.cpp

@@ -386 +386 @@
         /* chng 03 jan 13 rjrw, add in dynamics if required here,
          * last test - only do advection if we have not overrun the radius scale */
-        if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection 
+        if( dynamics.lgAdvection && iteration > dynamics.n_initial_relax+1 
                 /*&& radius.depth < dynamics.oldFullDepth*/ )
         {
@@ -758 +758 @@
 
         renormnew = 1.;
-        /* >>chng 06 mar 17, comment out test on old full depth - keep old solution if overrun scale */
-        if(iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection /*&& radius.depth < dynamics.oldFullDepth */)
-        {
-                /* The normalization out of the matrix solution is correct and
-                 * should be retained if: dynamics is on and the total
-                 * abundance of HI & H+ isn't being controlled by the
-                 * molecular network */
-                renormnew = 1.;
-        }
-        else
         {
                 double dennew = 0.,

branches/newmole/source/mole_newton_step.cpp

@@ -7 +7 @@
 /*lint -e725 expect positive indentation */
 #include "cddefines.h"
+#include "dense.h"
 #include "dynamics.h"
 #include "mole.h"
@@ -33 +34 @@
 
 static void updateMolecules(double *b2, double fion[LIMELM][N_MOLE_ION], realnum numelem[]);
-static void mole_eval_dynamic_balance(long int num_total, double *b, double **c);
+STATIC void mole_eval_dynamic_balance(long int num_total, double *b, double **c);
 STATIC int32 solve_system(double *a, double *b, long int n);
 
@@ -702 +703 @@
         return;
 }
-
-void mole_eval_dynamic_balance(long int num_total, double *b, double **c)
+STATIC void mole_eval_dynamic_balance(long int num_total, double *b, double **c)
 {
-        long int i;
+        long int i,
+                ion,
+                nelem;
+        double
+                source;
         
         DEBUG_ENTRY("mole_eval_dynamic_balance()");
+
         mole_eval_balance(num_total, b, c);
 
@@ -714 +719 @@
         {
                 /* Don't use conservation form in matrix solution -- dynamics rate terms make c[][] non-singular */             
+
                 for(i=0;i<mole.num_calc;i++) 
                 {
                         if (c)
                                 c[COmole[i]->index][COmole[i]->index] -= dynamics.Rate;
-                        b[COmole[i]->index] -= (COmole[i]->hevmol-dynamics.molecules[i])*dynamics.Rate;
+
+                        b[COmole[i]->index] -= COmole[i]->hevmol*dynamics.Rate;
+
+                        if (COmole[i]->n_nuclei != 1 || COmole[i]->nElec < 0)
+                        {
+                                b[COmole[i]->index] += dynamics.molecules[i]*dynamics.Rate;
+                        }
+                        else if (COmole[i]->nElec == 0)
+                        {
+                                for ( nelem=ipHYDROGEN; nelem < LIMELM; nelem++)
+                                        if (COmole[i]->nElem[nelem] == 1)
+                                                break;
+                                source = 0.0;
+                                for ( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion ) 
+                                {
+                                        source += dynamics.Source[nelem][ion];
+                                        if (ion >= N_MOLE_ION || element_list[nelem]->ipMl[ion] == -1)
+                                                source -= dense.xIonDense[nelem][ion]*dynamics.Rate;
+                                }
+                                b[COmole[i]->index] += source;
+                        }
                 }
         }