Changeset 2044 for branches/newmole

Timestamp:

05/11/08 14:54:02 (8 months ago)

Author:

rjrw

Message:

Move towards using Householder reflections to apply species
conservation constraints.

Presently optional coding switched off by "if (0)" uses projections to
confirm method detects conservation conservation constraints in error
vector.

Using the Householder vectors to project the away the conservation
nullspace from the matrix given to the linear solver should improve
stability.

Files:

• branches/newmole/source/mole_solve.cpp (modified) (2 diffs)

branches/newmole/source/mole_solve.cpp

@@ -18 +18 @@
 #include "h2.h"
 #include "newton_step.h"
+#include "householder.h"
 /* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element 
  * molecular network in Cloudy. Before this routine would predict negative abundances if 
@@ -441 +442 @@
         if (lgConserve)
         {
+                // Development code towards using Householder reflections to apply
+                // conservation constraints
                 if (0)
                 {
-                        double sum1;
+                        multi_arr<double,2> u(mole.num_calc,LIMELM);
+                        valarray<double> betah(LIMELM);
+                        for(i=0;i<mole.num_calc;i++)
+                        {
+                                for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                                {
+                                        u[i][nelem] = mole.list[i]->nElem[nelem];
+                                }
+                        }
                         for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
                         {
-                                sum1 = 0.;
+                                double beta;
+                                valarray<double> x(mole.num_calc-nelem),
+                                        y(mole.num_calc-nelem),
+                                        v(mole.num_calc-nelem), 
+                                        p(mole.num_calc-nelem);
+                                for(i=0;i<mole.num_calc-nelem;i++)
+                                {
+                                        x[i] = u[i+nelem][nelem];
+                                }
+                                householder(x,v,beta);
+                                betah[nelem] = beta;
+                                for(i=0;i<mole.num_calc-nelem;i++)
+                                {
+                                        u[i+nelem][nelem] = v[i];
+                                }
+                                for(long nother=nelem+1;nother<LIMELM;nother++)
+                                {                                
+                                        for (i=0;i<mole.num_calc-nelem;i++)
+                                        {
+                                                y[i] = u[i+nelem][nother];
+                                        }
+                                        p = project(y,v,beta);
+                                        for (i=0;i<mole.num_calc-nelem;i++)
+                                        {
+                                                u[i+nelem][nother] = p[i];
+                                        }                                       
+                                }
+                        }
+                        valarray<double> b2(mole.num_calc);
+                        for(i=0;i<mole.num_calc;i++)
+                        {
+                                b2[i] = b[i];
+                        }
+                        for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                        {
+                                valarray<double> v(mole.num_calc-nelem),
+                                        c(mole.num_calc-nelem),
+                                        p(mole.num_calc-nelem);
+                                for(i=0;i<mole.num_calc-nelem;i++)
+                                {
+                                        v[i] = u[i+nelem][nelem];
+                                        c[i] = b2[i+nelem];
+                                }
+                                p = project(c,v,betah[nelem]);
+                                for (i=0;i<mole.num_calc-nelem;i++)
+                                {
+                                        b2[i+nelem] = p[i];
+                                }
+                        }
+                        for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                        {
+                                fprintf(ioQQQ,"Elem %ld proj %g %g\n",nelem,b2[nelem],dense.gas_phase[nelem]);
+                        }
+                        for (;nelem<mole.num_calc;nelem++)
+                        {
+                                fprintf(ioQQQ,"Remain %ld proj %g\n",nelem,b2[nelem]);
+                        }
+                }
+                if (0)
+                {
+                        for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
+                        {
+                                double sum1 = 0.;
                                 for(i=0;i<mole.num_calc;i++)
                                 {