Context Navigation

← Previous Change
Next Change →

source

Timestamp:

02/02/08 10:58:59 (10 months ago)

Author:

rjrw

Message:

Complete factoring chemical solver detail out of mole_newton_step,
rename to newton_step.

Still cares about trying to get positive solution.

Location:

branches/newmole/source

Files:

: 1 added
: 3 modified
: 1 moved

cont_createmesh.cpp (modified) (1 diff)
mole_priv.h (modified) (2 diffs)
mole_solve.cpp (modified) (4 diffs)
newton_step.cpp (moved) (moved from branches/newmole/source/mole_newton_step.cpp) (16 diffs)
newton_step.h (added)

Legend:

: Unmodified
: Added
: Removed

branches/newmole/source/cont_createmesh.cpp

r1739	r1769
646	646
647	647	strcpy( chFilename , cpu.DataPath() );
648		strcat( chFilename , "continuum_mesh.~~dat~~" );
	648	strcat( chFilename , "continuum_mesh.ini" );
649	649
650	650	if( trace.lgTrace )

branches/newmole/source/mole_priv.h

r1768	r1769
13	13	map <string,double ()(struct mole_reaction_s rate)> functab;
14	14	} mole_priv;
	15
	16	class GroupMap {
	17	public:
	18	double fion[LIMELM][N_MOLE_ION];
	19	void updateMolecules(double *b2);
	20	void setup(double *b0vec);
	21	};
	22
	23	extern GroupMap MoleMap;
	24
15	25	typedef map<string,struct mole_reaction_s *>::iterator mole_reaction_i;
16	26	typedef map<string,struct molecule *>::iterator molecule_i;
…	…
58	68	extern void mole_update_rks( void );
59	69
60		/** Take one Newton step ~~of the chemical network~~
	70	/** Take one Newton step
61	71	\param *nBad
62	72	\param *error
63	73	*/
64		void ~~mole_newton_step(~~int nBad, realnum error, long n);
	74	void newton_step(valarray<double> b2vec, int nBad, realnum error, long n);
65	75
66	76	#endif /* _MOLE_PRIV_H_ */

branches/newmole/source/mole_solve.cpp

r1768	r1769
17	17	#include "grainvar.h"
18	18	#include "h2.h"
	19	#include "newton_step.h"
19	20	/* Nick Abel between July and October of 2003 assisted Dr. Ferland in improving the heavy element
20	21	* molecular network in Cloudy. Before this routine would predict negative abundances if
…	…
42	43	realnum
43	44	error;
	45	valarray<double> b2vec(mole.num_compacted);
44	46
45	47	DEBUG_ENTRY( "mole_solve()" );
…	…
88	90	nPrevBad = nBad;
89	91
90		mole_newton_step(&nBad, &error, mole.num_compacted);
	92	MoleMap.setup(&b2vec[0]);
	93	newton_step(b2vec, &nBad, &error, mole.num_compacted);
	94	MoleMap.updateMolecules( &b2vec[0] );
91	95
92	96	//fprintf(stdout,"Mole zone %ld -- %7.2f loop %d error %15.8g nBad %d\n",nzone,fnzone,i,error,nBad);
…	…
264	268	}
265	269
	270	#define ABSLIM 1e-12
	271	#define ERRLIM 1e-12
	272	#define SMALLABUND 1e-24
	273	# ifdef MAT
	274	# undef MAT
	275	# endif
	276	# define MAT(a,I_,J_) ((a)[(I_)*(mole.num_compacted)+(J_)])
	277
	278
	279	STATIC void mole_eval_dynamic_balance(long int num_total, double b, double *c);
	280	enum {PRINTSOL = false};
	281
	282	GroupMap MoleMap;
	283
	284	void funjac(valarray<double> b2vec, double ervals, double amat, bool lgJac)
	285	{
	286	long
	287	nelem,
	288	ion,
	289	i,
	290	j;
	291	double sum;
	292	static double
	293	**c;
	294	int printsol = PRINTSOL;
	295	valarray<double> b(mole.num_total);
	296	static bool lgMustMalloc = true;
	297	bool lgConserve;
	298
	299	DEBUG_ENTRY( "funjac()" );
	300
	301	MoleMap.updateMolecules( &b2vec[0] );
	302
	303	if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection /&& radius.depth < dynamics.oldFullDepth/ )
	304	{
	305	ASSERT(dynamics.Rate > 0.0);
	306	lgConserve = false;
	307	}
	308	else
	309	{
	310	lgConserve = true;
	311	}
	312
	313	if( lgMustMalloc )
	314	{
	315	/* on very first call must create space */
	316	lgMustMalloc = false;
	317	c = (double *)MALLOC((size_t)mole.num_totalsizeof(double *));
	318	c[0] = (double )MALLOC((size_t)mole.num_total
	319	mole.num_total*sizeof(double));
	320	for(i=1;i<mole.num_total;i++)
	321	{
	322	c[i] = c[i-1]+mole.num_total;
	323	}
	324	}
	325
	326	/* Generate chemical balance vector (mole.b[]) and Jacobian array
	327	(c[][], first iteration) from reaction list */
	328
	329	mole_eval_dynamic_balance(mole.num_total,&b[0],lgJac?c:NULL);
	330
	331	/* add positive ions and neutral atoms: ratios are set by ion_solver,
	332	* we determine abundance of the group as a whole here */
	333
	334	if (lgJac) {
	335	// c[i][i] can be positive for auto-catalytic reactions
	336	//for(i=0;i<mole.num_calc;i++)
	337	//{
	338	// ASSERT (c[i][i] <= 0.);
	339	//}
	340	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	341	{
	342	if (element_list[nelem]->ipMl[0] != -1)
	343	{
	344	for(i=0;i<mole.num_calc;i++)
	345	{
	346	ASSERT(mole.list[i]->index == i);
	347	sum = 0.;
	348	for (ion=0;ion<N_MOLE_ION;ion++)
	349	{
	350	if (element_list[nelem]->ipMl[ion] != -1)
	351	{
	352	sum += MoleMap.fion[nelem][ion]*c[element_list[nelem]->ipMl[ion]][i];
	353	c[element_list[nelem]->ipMl[ion]][i] = 0.;
	354	}
	355	}
	356	c[element_list[nelem]->ipMl[0]][i] = sum;
	357	}
	358	}
	359	}
	360	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	361	{
	362	if (element_list[nelem]->ipMl[0] != -1)
	363	{
	364	for(i=0;i<mole.num_calc;i++)
	365	{
	366	sum = 0.0;
	367	for (ion=0;ion<N_MOLE_ION;ion++)
	368	{
	369	if (element_list[nelem]->ipMl[ion] != -1)
	370	{
	371	sum += c[i][element_list[nelem]->ipMl[ion]];
	372	c[i][element_list[nelem]->ipMl[ion]] = 0.;
	373	}
	374	}
	375	c[i][element_list[nelem]->ipMl[0]] = sum;
	376	}
	377	}
	378	}
	379	}
	380
	381	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	382	{
	383	if (element_list[nelem]->ipMl[0] != -1)
	384	{
	385	sum = 0.0;
	386	for (ion=0;ion<N_MOLE_ION;ion++)
	387	{
	388	if (element_list[nelem]->ipMl[ion] != -1)
	389	{
	390	sum += b[element_list[nelem]->ipMl[ion]];
	391	b[element_list[nelem]->ipMl[ion]] = 0.0;
	392	}
	393	}
	394	b[element_list[nelem]->ipMl[0]] = sum;
	395	}
	396	}
	397
	398	/* For Newton-Raphson method, want the change in populations to be zero,
	399	* so the conserved component must also be zero */
	400	if (lgConserve)
	401	{
	402	double scale;
	403	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	404	{
	405	if (element_list[nelem]->ipMl[0] != -1)
	406	{
	407	if (lgJac)
	408	{
	409	scale = -(ABSLIM+fabs(c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]]));
	410	for(i=0;i<mole.num_calc;i++)
	411	{
	412	c[i][element_list[nelem]->ipMl[0]] = mole.list[i]->nElem[nelem]*scale;
	413	}
	414	}
	415	else
	416	scale = c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]];
	417
	418	b[element_list[nelem]->ipMl[0]] = 0.;
	419	}
	420	}
	421	}
	422
	423	/------------------------------------------------------------------ /
	424	if(printsol \|\| (trace.lgTrace && trace.lgTr_H2_Mole ))
	425	{
	426	/* print the full matrix */
	427	fprintf( ioQQQ, " ");
	428	for( i=0; i < mole.num_calc; i++ )
	429	{
	430	fprintf( ioQQQ, " %-4.4s", mole.list[i]->label );
	431	}
	432	fprintf( ioQQQ, " \n" );
	433
	434	fprintf(ioQQQ," MOLE old abundances\t%.2f",fnzone);
	435	for( i=0; i<mole.num_calc; i++ )
	436	fprintf(ioQQQ,"\t%.2e", mole.list[i]->den );
	437	fprintf(ioQQQ,"\n" );
	438
	439	for( i=0; i < mole.num_calc; i++ )
	440	{
	441	fprintf( ioQQQ, " MOLE%2ld %-4.4s", i ,mole.list[i]->label);
	442	for( j=0; j < mole.num_calc; j++ )
	443	{
	444	fprintf( ioQQQ, "%10.2e", c[j][i] );
	445	}
	446	fprintf( ioQQQ, "%10.2e", b[i] );
	447	fprintf( ioQQQ, "\n" );
	448	}
	449	}
	450
	451	for( i=0; i < mole.num_compacted; i++ )
	452	{
	453	ervals[i] = b[groupspecies[i]->index];
	454	}
	455
	456	if (lgJac)
	457	{
	458	for( i=0; i < mole.num_compacted; i++ )
	459	{
	460	// c[i][i] can be +ve for (effectively) autocatalytic reactions,
	461	// e.g. H-,H+=>H,H
	462	//ASSERT (c[groupspecies[i]->index][groupspecies[i]->index] <= 0.);
	463	for( j=0; j < mole.num_compacted; j++ )
	464	{
	465	MAT(amat,i,j) = c[groupspecies[i]->index][groupspecies[j]->index];
	466	}
	467	}
	468	}
	469	}
	470
	471	void GroupMap::setup(double *b0vec)
	472	{
	473	valarray<double> calcv(mole.num_total);
	474	long i,
	475	nelem,
	476	ion;
	477	bool lgSet;
	478	double sum;
	479
	480	for(i=0;i<mole.num_total;i++)
	481	{
	482	calcv[i] = mole.list[i]->den;
	483	}
	484
	485	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	486	{
	487	if (element_list[nelem]->ipMl[0] != -1)
	488	{
	489	sum = 0.;
	490	for (ion=0; ion<N_MOLE_ION; ion++)
	491	{
	492	if (element_list[nelem]->ipMl[ion] != -1)
	493	sum += calcv[element_list[nelem]->ipMl[ion]];
	494	}
	495	if (sum > SMALLFLOAT)
	496	{
	497	for (ion=0; ion<N_MOLE_ION; ion++)
	498	{
	499	if (element_list[nelem]->ipMl[ion] != -1)
	500	fion[nelem][ion] = calcv[element_list[nelem]->ipMl[ion]]/sum;
	501	else
	502	fion[nelem][ion] = 0.;
	503	}
	504	}
	505	else
	506	{
	507	lgSet = false;
	508	for (ion=0; ion<N_MOLE_ION; ion++)
	509	{
	510	if (element_list[nelem]->ipMl[ion] != -1 && !lgSet)
	511	{
	512	fion[nelem][ion] = 1.0;
	513	lgSet = true;
	514	}
	515	else
	516	{
	517	fion[nelem][ion] = 0.;
	518	}
	519	}
	520	}
	521
	522	lgSet = false;
	523	for (ion=0;ion<N_MOLE_ION;ion++)
	524	{
	525	if (element_list[nelem]->ipMl[ion] != -1)
	526	{
	527	if (!lgSet)
	528	calcv[element_list[nelem]->ipMl[ion]] = sum;
	529	else
	530	calcv[element_list[nelem]->ipMl[ion]] = 0.;
	531	lgSet = true;
	532	}
	533	}
	534	}
	535	}
	536
	537	for( i=0; i < mole.num_compacted; i++ )
	538	{
	539	b0vec[i] = calcv[groupspecies[i]->index];
	540	}
	541	}
	542
	543	void GroupMap::updateMolecules(double *b2)
	544	{
	545	long int
	546	ion,
	547	mol,
	548	nelem;
	549	double
	550	grouptot,
	551	sum;
	552
	553	DEBUG_ENTRY("updateMolecules()");
	554
	555	for (mol=0;mol<mole.num_calc;mol++)
	556	{
	557	mole.list[mol]->den = 0.;
	558	}
	559	for (mol=0;mol<mole.num_compacted;mol++)
	560	{
	561	groupspecies[mol]->den = b2[mol]; /* put derived abundances back into appropriate molecular species */
	562	}
	563
	564	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	565	{
	566	if (element_list[nelem]->ipMl[0] != -1)
	567	{
	568	grouptot = mole.list[element_list[nelem]->ipMl[0]]->den;
	569	sum = 0.0;
	570	for (ion=0;ion<N_MOLE_ION;ion++)
	571	{
	572	if (element_list[nelem]->ipMl[ion] != -1)
	573	{
	574	mole.list[element_list[nelem]->ipMl[ion]]->den = grouptot * fion[nelem][ion];
	575	sum += mole.list[element_list[nelem]->ipMl[ion]]->den;
	576	}
	577	}
	578	ASSERT(fabs(sum-grouptot) <= 1e-10 * fabs(grouptot));
	579	}
	580	}
	581
	582	return;
	583	}
	584
	585	STATIC void mole_eval_dynamic_balance(long int num_total, double b, double *c)
	586	{
	587	long int i,
	588	ion,
	589	nelem;
	590	double
	591	source;
	592
	593	DEBUG_ENTRY("mole_eval_dynamic_balance()");
	594
	595	mole_eval_balance(num_total, b, c);
	596
	597	/* >>chng 06 mar 17, comment out test on old full depth - keep old solution if overrun scale */
	598	if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection /&& radius.depth < dynamics.oldFullDepth/ )
	599	{
	600	/* Don't use conservation form in matrix solution -- dynamics rate terms make c[][] non-singular */
	601
	602	for(i=0;i<mole.num_calc;i++)
	603	{
	604	if (c)
	605	c[i][i] -= dynamics.Rate;
	606
	607	b[i] -= mole.list[i]->den*dynamics.Rate;
	608
	609	if (mole.list[i]->n_nuclei != 1 \|\| mole.list[i]->charge < 0 \|\| ! mole.list[i]->lgGas_Phase)
	610	{
	611	b[i] += dynamics.molecules[i]*dynamics.Rate;
	612	}
	613	else if (mole.list[i]->charge == 0)
	614	{
	615	for ( nelem=ipHYDROGEN; nelem < LIMELM; nelem++)
	616	if (mole.list[i]->nElem[nelem] == 1)
	617	break;
	618	source = 0.0;
	619	for ( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
	620	{
	621	source += dynamics.Source[nelem][ion];
	622	if (ion >= N_MOLE_ION \|\| element_list[nelem]->ipMl[ion] == -1)
	623	source -= dense.xIonDense[nelem][ion]*dynamics.Rate;
	624	}
	625	b[i] += source;
	626	}
	627	}
	628	}
	629	}

branches/newmole/source/newton_step.cpp

r1768	r1769
7	7	/lint -e725 expect positive indentation /
8	8	#include "cddefines.h"
9		~~#include "dense.h"~~
10		~~#include "dynamics.h"~~
11		~~#include "mole.h"~~
12		~~#include "mole_priv.h"~~
13		~~#include "physconst.h"~~
14	9	#include "thirdparty.h"
15		~~#include "trace.h"~~
16
17		~~/* >>chng 07 dec 10 rjrw -- it's simpler now */~~
18		~~/* HP cc cannot compile following except in -O1 mode */~~
19		~~//#if defined(__HP_aCC)~~
20		~~//# pragma OPT_LEVEL 1~~
21		~~//#endif~~
22		~~// icc 10.0 miscompiles this routine with higher optimization~~
23		~~//#if defined (__ICC) && defined(__amd64)~~
24		~~//#pragma optimization_level 1~~
25		~~//#endif~~
26	10
27	11	#define ABSLIM 1e-12
…	…
31	15	# undef MAT
32	16	# endif
33		# define MAT(a,I_,J_) ((a)[(I_)*(mole.num_compacted)+(J_)])
34
35
36		STATIC void mole_eval_dynamic_balance(long int num_total, double b, double *c);
37		STATIC int32 solve_system(double a, double b, long int n);
	17	# define MAT(a,I_,J_) ((a)[(I_)*(n)+(J_)])
	18
	19
	20	STATIC int32 solve_system(const valarray<double> &a, valarray<double> &b, const long int n);
38	21	enum {PRINTSOL = false};
39	22
40		static class GroupMap {
41		public:
42		double fion[LIMELM][N_MOLE_ION];
43		void updateMolecules(double *b2);
44		void setup(double *b0vec);
45		} MoleMap;
46
47		static void funjac(double ervals, double amat, long loop)
48		{
49		long
50		nelem,
51		ion,
52		i,
53		j;
54		double sum;
55		static double
56		**c;
57		int printsol = PRINTSOL;
58		valarray<double> b(mole.num_total);
59		static bool lgMustMalloc = true;
60		bool lgConserve;
61
62		DEBUG_ENTRY( "funjac()" );
63
64		if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection /&& radius.depth < dynamics.oldFullDepth/ )
65		{
66		ASSERT(dynamics.Rate > 0.0);
67		lgConserve = false;
68		}
69		else
70		{
71		lgConserve = true;
72		}
73
74		if( lgMustMalloc )
75		{
76		/* on very first call must create space */
77		lgMustMalloc = false;
78		c = (double *)MALLOC((size_t)mole.num_totalsizeof(double *));
79		c[0] = (double )MALLOC((size_t)mole.num_total
80		mole.num_total*sizeof(double));
81		for(i=1;i<mole.num_total;i++)
82		{
83		c[i] = c[i-1]+mole.num_total;
84		}
85		}
86
87		/* Generate chemical balance vector (mole.b[]) and Jacobian array
88		(c[][], first iteration) from reaction list */
89
90		mole_eval_dynamic_balance(mole.num_total,&b[0],(loop==0)?c:NULL);
91
92		/* add positive ions and neutral atoms: ratios are set by ion_solver,
93		* we determine abundance of the group as a whole here */
94
95		if (loop == 0) {
96		// c[i][i] can be positive for auto-catalytic reactions
97		//for(i=0;i<mole.num_calc;i++)
98		//{
99		// ASSERT (c[i][i] <= 0.);
100		//}
101		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
102		{
103		if (element_list[nelem]->ipMl[0] != -1)
104		{
105		for(i=0;i<mole.num_calc;i++)
106		{
107		ASSERT(mole.list[i]->index == i);
108		sum = 0.;
109		for (ion=0;ion<N_MOLE_ION;ion++)
110		{
111		if (element_list[nelem]->ipMl[ion] != -1)
112		{
113		sum += MoleMap.fion[nelem][ion]*c[element_list[nelem]->ipMl[ion]][i];
114		c[element_list[nelem]->ipMl[ion]][i] = 0.;
115		}
116		}
117		c[element_list[nelem]->ipMl[0]][i] = sum;
118		}
119		}
120		}
121		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
122		{
123		if (element_list[nelem]->ipMl[0] != -1)
124		{
125		for(i=0;i<mole.num_calc;i++)
126		{
127		sum = 0.0;
128		for (ion=0;ion<N_MOLE_ION;ion++)
129		{
130		if (element_list[nelem]->ipMl[ion] != -1)
131		{
132		sum += c[i][element_list[nelem]->ipMl[ion]];
133		c[i][element_list[nelem]->ipMl[ion]] = 0.;
134		}
135		}
136		c[i][element_list[nelem]->ipMl[0]] = sum;
137		}
138		}
139		}
140		}
141
142		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
143		{
144		if (element_list[nelem]->ipMl[0] != -1)
145		{
146		sum = 0.0;
147		for (ion=0;ion<N_MOLE_ION;ion++)
148		{
149		if (element_list[nelem]->ipMl[ion] != -1)
150		{
151		sum += b[element_list[nelem]->ipMl[ion]];
152		b[element_list[nelem]->ipMl[ion]] = 0.0;
153		}
154		}
155		b[element_list[nelem]->ipMl[0]] = sum;
156		}
157		}
158
159		/* For Newton-Raphson method, want the change in populations to be zero,
160		* so the conserved component must also be zero */
161		if (lgConserve)
162		{
163		double scale;
164		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
165		{
166		if (element_list[nelem]->ipMl[0] != -1)
167		{
168		if (loop == 0)
169		{
170		scale = -(ABSLIM+fabs(c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]]));
171		for(i=0;i<mole.num_calc;i++)
172		{
173		c[i][element_list[nelem]->ipMl[0]] = mole.list[i]->nElem[nelem]*scale;
174		}
175