Changeset 2000

Timestamp:

04/30/08 16:02:48 (2 weeks ago)

Author:

rjrw

Message:

Initial read across of ion sources to isolevel sources.

Runs the smoke test without crashing, now need to correct for the
different basis of StatesElem? compared to xIonDense.

Files:

• branches/isodyn/source/atmdat_readin.cpp (modified) (2 diffs)
• branches/isodyn/source/dynamics.cpp (modified) (15 diffs)
• branches/isodyn/source/radius_increment.cpp (modified) (1 diff)
• branches/isodyn/source/state.cpp (modified) (1 diff)
• branches/isodyn/source/struc.h (modified) (1 diff)

branches/isodyn/source/atmdat_readin.cpp

@@ -21 +21 @@
-#include "hyperfine.h"
-#include "atmdat.h"
+#include "iso.h"
-/* */
-/* this was needed to get array to crash out of bounds if not set.
@@ -303 +304 @@
-        }
-
+        struc.StatesElem = (realnum ****)MALLOC(sizeof(realnum ***)*(unsigned)(NISO) );
+        for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
+        {
+                struc.StatesElem[ipISO] = 
+                        (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(LIMELM) );
+                for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
+                {
+                        struc.StatesElem[ipISO][nelem] = 
+                                (realnum**)MALLOC(sizeof(realnum*)*(unsigned)iso.numLevels_max[ipISO][nelem]);
+                        for( long level=0; level < iso.numLevels_max[ipISO][nelem]; ++level )
+                        {
+                                struc.StatesElem[ipISO][nelem][level] = 
+                                        (realnum*)MALLOC(sizeof(realnum)*(unsigned)(struc.nzlim) );
+                        }
+                }
+        }
+
-        /* some structure variables */
-        for( i=0; i < struc.nzlim; i++ )

branches/isodyn/source/dynamics.cpp

@@ -36 +36 @@
-#include "wind.h"
-#include "hmi.h"
+#include "iso.h"
-#include "dynamics.h"
-static int ipUpstream=0,iphUpstream=0,ipyUpstream=0;
@@ -57 +58 @@
- * [element][ion] */
-static double **UpstreamIon;
+static double ***UpstreamStatesElem;
-/* total abundance of each element per hydrogen */
-static double *UpstreamElem;
@@ -135 +137 @@
-/* the gas phase abundances from the previous iteration */
-static realnum **Old_gas_phase;
+
+/* the iso levels from the previous iteration */
+static realnum ****Old_StatesElem;
-
-/* the number of zones that were saved in the previous iteration */
@@ -840 +845 @@
-
-        /* Properties for cell half as far upstream, used to converge timestep */
-
+, hiso
-
-        /* Properties for cell at present depth, used to converge timestep */
-
+, yiso
-
-        long int nelem , ion, mol;
@@ -859 +864 @@
-                        {
-                                UpstreamIon[nelem][ion] = 0.;
+                        }
+                }
+                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 )
+                                {
+                                        UpstreamStatesElem[ipISO][nelem][level] = 0.;
+                                }
-                        }
-                }
@@ -955 +970 @@
-                }
-
+                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 )
+                                {
+                                        UpstreamStatesElem[ipISO][nelem][level] = 
+                                                Old_StatesElem[ipUpstream][ipISO][nelem][level]/Old_hden[ipUpstream] +
+                                                (Old_StatesElem[ipUpstream+1][ipISO][nelem][level]/Old_hden[ipUpstream+1] - 
+                                                 Old_StatesElem[ipUpstream][ipISO][nelem][level]/Old_hden[ipUpstream])*
+                                                frac_next;
+                                }
+                        }
+                }
+
-                for(mol=0;mol<N_H_MOLEC;mol++)
-                {
@@ -1005 +1035 @@
-                                        Old_xIonDense[ipUpstream][nelem][ion]/Old_hden[ipUpstream];
-                                UpstreamElem[nelem] += UpstreamIon[nelem][ion];
+                        }
+                }
+
+                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 )
+                                {
+                                        UpstreamStatesElem[ipISO][nelem][level] = 
+                                                Old_StatesElem[ipUpstream][ipISO][nelem][level]/Old_hden[ipUpstream];
+                                }
-                        }
-                }
@@ -1134 +1176 @@
-                }
-        }
+
+        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 )
+                        {
+                                double f1;
+                                if(ipyUpstream != nOld_zone-1 && (Old_depth[ipyUpstream+1] - Old_depth[ipyUpstream]>SMALLFLOAT))
+                                {
+                                        yiso = 
+                                                Old_StatesElem[ipyUpstream][ipISO][nelem][level]/Old_hden[ipyUpstream] +
+                                                (Old_StatesElem[ipyUpstream+1][ipISO][nelem][level]/Old_hden[ipyUpstream+1] - 
+                                                 Old_StatesElem[ipyUpstream][ipISO][nelem][level]/Old_hden[ipyUpstream])*
+                                                ynextfrac;
+                                }
+                                else
+                                {               
+                                        yiso = Old_StatesElem[ipyUpstream][ipISO][nelem][level]/Old_hden[ipyUpstream];
+                                }
+                                if(iphUpstream != nOld_zone-1 && (Old_depth[iphUpstream+1] - Old_depth[iphUpstream]>SMALLFLOAT))
+                                {
+                                        hiso = 
+                                                Old_StatesElem[iphUpstream][ipISO][nelem][level]/Old_hden[iphUpstream] +
+                                                (Old_StatesElem[iphUpstream+1][ipISO][nelem][level]/Old_hden[iphUpstream+1] - 
+                                                 Old_StatesElem[iphUpstream][ipISO][nelem][level]/Old_hden[iphUpstream])*
+                                                hnextfrac;
+                                }
+                                else
+                                {               
+                                        hiso = Old_StatesElem[iphUpstream][ipISO][nelem][level]/Old_hden[iphUpstream];
+                                }
+
+                                /* the proposed thickness of the next zone, there will be a scale factor, something
+                                 * like 1/500, that will be applied when this is used in nextdr */
+                                f1 = fabs(yiso - UpstreamStatesElem[ipISO][nelem][level] );
+                                f1 = SDIV( f1 );
+                                dynamics.dRad = MIN2(dynamics.dRad,Dyn_dr * 
+                          /* don't pay attention to species with abundance relative to H less tghan 1e-6 */
+                                  MAX2(yiso + UpstreamStatesElem[ipISO][nelem][level],1e-6 ) / f1);
+                                /*MAX2(fabs(yion - UpstreamIon[nelem][ion] ),SMALLFLOAT));*/
+                                /* printf("%d %d %g\n",nelem,ion,dynamics.dRad); */
+
+                                /* Must be consistent with convergence_error below */
+                                /* >>chngf 01 aug 01, remove dense.gas_phase[ipHYDROGEN] from HAtom */
+                                /* >>chngf 02 aug 01, multiply by cell width */
+                                /* this error is error due to the advection length not being zero - a finite
+                                 * advection length.  no need to bring convergence error to below
+                                 * discretization error.  when convergece error is lower than a fraction of this
+                                 * errror we reduce the advection length. */
+                                dynamics.discretization_error += POW2(yiso-2.*hiso+UpstreamStatesElem[ipISO][nelem][level]); 
+                                dynamics.error_scale2 += POW2(UpstreamStatesElem[ipISO][nelem][level]);
+                        }
+                }
+        }
+
-
-        for(mol=0; mol < N_H_MOLEC; mol++) 
@@ -1542 +1640 @@
-                Oldi_hden,
-                Oldi_ion,
+                Oldi_iso,
-                Oldi_mol;
-
@@ -1593 +1692 @@
-                                /* >>chng 02 nov 11, add first error scale estimate from Robin */
-                                dynamics.error_scale1 += POW2(struc.xIonDense[nelem][ion][i]/struc.hden[i]);                    
+                        }
+                }
+                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 )
+                                {
+                                        if(ilast != nOld_zone-1 && ((Old_depth[ilast+1] - Old_depth[ilast])> SMALLFLOAT) )
+                                        {
+                                                Oldi_iso = (Old_StatesElem[ilast][ipISO][nelem][level] +
+                                                  (Old_StatesElem[ilast+1][ipISO][nelem][level]-Old_StatesElem[ilast][ipISO][nelem][level])*
+                                                                                                frac_next);
+                                        }
+                                        else
+                                        {
+                                                Oldi_iso = Old_StatesElem[ilast][ipISO][nelem][level];
+                                        }
+                                        dynamics.convergence_error += POW2(Oldi_iso/Oldi_hden-struc.StatesElem[ipISO][nelem][level][i]/struc.hden[i]) /* *struc.dr[i] */;  
+
+                                        /* >>chng 02 nov 11, add first error scale estimate from Robin */
+                                        dynamics.error_scale1 += POW2(struc.StatesElem[ipISO][nelem][level][i]/struc.hden[i]);                  
+                                }
-                        }
-                }
@@ -1698 +1820 @@
-                        }
-                }
+                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 )
+                                { 
+                                        Old_StatesElem[i][ipISO][nelem][level] = struc.StatesElem[ipISO][nelem][level][i];
+                                }
+                        }
+                }               
-        }
-        return;
@@ -1861 +1993 @@
-        Old_xIonDense = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)(struc.nzlim) );
-
+        Old_StatesElem = (realnum ****)MALLOC(sizeof(realnum ***)*(unsigned)(struc.nzlim) );
+
-        Old_gas_phase = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)(struc.nzlim) );
-
@@ -1873 +2007 @@
-                        (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM+3) );
-
+                Old_StatesElem[ns] = 
+                        (realnum***)MALLOC(sizeof(realnum**)*(unsigned)(NISO) );
+
-                Old_gas_phase[ns] = 
-                        (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+3) );
@@ -1886 +2023 @@
-                        Old_xIonDense[ns][nelem] = 
-                                (realnum*)MALLOC(sizeof(realnum)*(unsigned)(LIMELM+1) );
+                }
+                for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
+                {
+                        Old_StatesElem[ns][ipISO] = 
+                                (realnum**)MALLOC(sizeof(realnum*)*(unsigned)(LIMELM) );
+                        for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem)
+                        {
+                                Old_StatesElem[ns][ipISO][nelem] = 
+                                        (realnum*)MALLOC(sizeof(realnum)*(unsigned)iso.numLevels_max[ipISO][nelem]);
+                        }
-                }
-        }
@@ -1910 +2057 @@
-                        {
-                                Old_xIonDense[i][nelem][ion] = 0.;
+                        }
+                }
+                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 )
+                                {
+                                        Old_StatesElem[i][ipISO][nelem][level] = 0.;
+                                }
-                        }
-                }

branches/isodyn/source/radius_increment.cpp

@@ -727 +727 @@
-                }
-        }
+        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 )
+                        {
+                                struc.StatesElem[ipISO][nelem][level][nzone_minus_1] = StatesElem[ipISO][nelem][level].Pop;
+                        }
+                }
+        }
-
-        /* the hydrogen molecules */

branches/isodyn/source/state.cpp

@@ -324 +324 @@
-                }
-        }
+        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 )
+                        {
+                                state_do( struc.StatesElem[ipISO][nelem][level] , (size_t)(struc.nzlim)*sizeof(realnum ) );
+                        }
+                }
+        }
-
-        /* the hydrogen molecules */

branches/isodyn/source/struc.h

@@ -62 +62 @@
-        realnum ***xIonDense;
-
+        /** save iso level array across model */
+        realnum ****StatesElem;
+
-        /** the hydrogen molecules */
-        /**realnum *Molec[N_H_MOLEC];*/