Changeset 2011

Timestamp:

05/02/08 16:50:16 (1 week ago)

Author:

rjrw

Message:

First cut at implementing dynamics terms.

Passes smoke test, fails during allocation for dynamics run (trying to
allocate before the size of the isolevel arrays is determined).

Files:

• branches/isodyn/source/dynamics.cpp (modified) (4 diffs)
• branches/isodyn/source/dynamics.h (modified) (1 diff)
• branches/isodyn/source/iso_level.cpp (modified) (2 diffs)
• branches/isodyn/source/service.cpp (modified) (1 diff)

branches/isodyn/source/dynamics.cpp

@@ -652 +652 @@
-                        }
-                }
-
+
+
+
+
+
+
+
+
+
+
-                for(mol=0;mol<N_H_MOLEC;mol++)
-                {
@@ -791 +800 @@
-                }
-        }
+
+        for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
+        {
+                for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
+                {
+                        for( long level=0; level < iso.numLevels_local[ipISO][nelem]; ++level )
+                        {
+                                dynamics.StatesElem[ipISO][nelem][level] = 
+                                        UpstreamStatesElem[ipISO][nelem][level] * dense.gas_phase[ipHYDROGEN] / timestep;
+                        }
+                }
+        }
+
-#       if 0
-        fprintf(ioQQQ,"dynamiccc\t%li\t%.2e\t%.2e\t%.2e\t%.2e\n",
@@ -1826 +1848 @@
-                                for( long level=0; level < iso.numLevels_local[ipISO][nelem]; ++level )
-                                { 
-
+*struc.xIonDense[nelem][nelem-ipISO][i]
-                                }
-                        }
@@ -1969 +1991 @@
-                }
-        }
+
+        dynamics.StatesElem = ((double***)MALLOC( (size_t)NISO*sizeof(double **) ));
+
+        for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
+        {
+                dynamics.StatesElem[ipISO] = 
+                        (double**)MALLOC(sizeof(double*)*(unsigned)(LIMELM) );
+                for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
+                {
+                        dynamics.StatesElem[ipISO][nelem] = 
+                                (double*)MALLOC(sizeof(double)*(unsigned)iso.numLevels_max[ipISO][nelem]);
+                }
+        }
+
-        dynamics.Rate = 0.;
-

branches/isodyn/source/dynamics.h

@@ -70 +70 @@
-        double Rate;
-
-recombination rate
+ionization balance terms
-        double **Source /*[LIMELM][LIMELM+2]*/;
+
+        /** the advective isolevel balance terms */
+        double ***StatesElem;
-
-        /** save H2 and CO densities */

branches/isodyn/source/iso_level.cpp

@@ -321 +321 @@
-                }
-
-                /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
-                /*>>>chng 02 oct 01, only include if lgAdvection is set */
-                if( dynamics.lgAdvection && dynamics.lgISO[ipISO])
-                {
-                        /* add in advection - these terms normally zero */
-                        /* assume for present that all advection is into ground state */
-                        source += dynamics.Source[nelem][nelem-ipISO];
-                        /* >>chng 02 Sep 06 rjrw -- advective term not recombination */
-                        /* can sink from all components (must do, for conservation) */
-                        sink += dynamics.Rate;
-                }
-
-#if     0
-                /* add in source and sink terms from molecular network. */
@@ -349 +337 @@
-                {
-                        z[level][level] += sink;
+                }
+
+                /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
+                /*>>>chng 02 oct 01, only include if lgAdvection is set */
+                if( dynamics.lgAdvection && dynamics.lgISO[ipISO])
+                {
+                        for( level=0; level < numlevels_local; level++ )
+                        {
+                                creation[level] += dynamics.StatesElem[nelem][ipISO][level]/SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
+                                z[level][level] += dynamics.Rate;
+                        }
-                }
-

branches/isodyn/source/service.cpp

@@ -1580 +1580 @@
-        DEBUG_ENTRY( "MyMalloc()" );
-
+        if (size == 0)
+                fprintf(stdout,"Failed at %s:%d\n",chFile,line);
-        ASSERT( size > 0 );
-