Context Navigation

← Previous Change
Next Change →

source

Timestamp:

02/01/08 15:14:18 (10 months ago)

Author:

rjrw

Message:

Start to factor chemistry out of Newton solver.

Location:

branches/newmole/source

Files:

: 3 modified

mole_newton_step.cpp (modified) (14 diffs)
mole_priv.h (modified) (1 diff)
mole_solve.cpp (modified) (1 diff)

Legend:

: Unmodified
: Added
: Removed

branches/newmole/source/mole_newton_step.cpp

r1741	r1768
31	31	# undef MAT
32	32	# endif
33		# define MAT(a,I_,J_) (((a)+(I_)(mole.num_compacted)+(J_)))
34
35
36		STATIC void updateMolecules(double *b2, double fion[LIMELM][N_MOLE_ION]);
	33	# define MAT(a,I_,J_) ((a)[(I_)*(mole.num_compacted)+(J_)])
	34
	35
37	36	STATIC void mole_eval_dynamic_balance(long int num_total, double b, double *c);
38	37	STATIC int32 solve_system(double a, double b, long int n);
39
	38	enum {PRINTSOL = false};
	39
	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	}
	176	else
	177	scale = c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]];
	178
	179	b[element_list[nelem]->ipMl[0]] = 0.;
	180	}
	181	}
	182	}
	183
	184	/------------------------------------------------------------------ /
	185	if(printsol \|\| (trace.lgTrace && trace.lgTr_H2_Mole ))
	186	{
	187	/* print the full matrix */
	188	fprintf( ioQQQ, " ");
	189	for( i=0; i < mole.num_calc; i++ )
	190	{
	191	fprintf( ioQQQ, " %-4.4s", mole.list[i]->label );
	192	}
	193	fprintf( ioQQQ, " \n" );
	194
	195	fprintf(ioQQQ," MOLE old abundances\t%.2f",fnzone);
	196	for( i=0; i<mole.num_calc; i++ )
	197	fprintf(ioQQQ,"\t%.2e", mole.list[i]->den );
	198	fprintf(ioQQQ,"\n" );
	199
	200	for( i=0; i < mole.num_calc; i++ )
	201	{
	202	fprintf( ioQQQ, " MOLE%2ld %-4.4s", i ,mole.list[i]->label);
	203	for( j=0; j < mole.num_calc; j++ )
	204	{
	205	fprintf( ioQQQ, "%10.2e", c[j][i] );
	206	}
	207	fprintf( ioQQQ, "%10.2e", b[i] );
	208	fprintf( ioQQQ, "\n" );
	209	}
	210	}
	211
	212	for( i=0; i < mole.num_compacted; i++ )
	213	{
	214	ervals[i] = b[groupspecies[i]->index];
	215	}
	216
	217	if (loop == 0)
	218	{
	219	for( i=0; i < mole.num_compacted; i++ )
	220	{
	221	// c[i][i] can be +ve for (effectively) autocatalytic reactions,
	222	// e.g. H-,H+=>H,H
	223	//ASSERT (c[groupspecies[i]->index][groupspecies[i]->index] <= 0.);
	224	for( j=0; j < mole.num_compacted; j++ )
	225	{
	226	MAT(amat,i,j) = c[groupspecies[i]->index][groupspecies[j]->index];
	227	}
	228	}
	229	}
	230	}
40	231
41	232
42	233	/* mole_newton_step -- improve balance in chemical network along
43	234	* descent direction, step limited to ensure improvement */
44		void mole_newton_step(int nBad, realnum error)
	235	void mole_newton_step(int nBad, realnum error, long n)
45	236	{
46		~~enum {PRINTSOL = false};~~
47
48	237	long int i,
49		j,
50		~~nelem,~~
51		~~ion,~~
52	238	loop;
53
54		~~bool lgConserve,~~
55		~~lgSet;~~
56	239
57	240	int printsol = PRINTSOL;
…	…
63	246	emax,
64	247	fracneg,
65		~~sum,~~
66		~~fion[LIMELM][N_MOLE_ION],~~
67	248	fstep,
68	249	flast,
…	…
71	252	grad;
72	253
73		static double
74		*amat=NULL,
75		*b0vec=NULL,
76		*b1vec=NULL,
77		*b2vec=NULL,
78		*calcv=NULL,
79		*b=NULL,
80		**c=NULL;
	254	valarray<double>
	255	amat(n*n),
	256	b0vec(n),
	257	b1vec(n),
	258	b2vec(n),
	259	escale(n),
	260	ervals(n);
81	261	int32 merror;
82	262
83		~~/* if this is still true then must create space for arrays */~~
84		~~static bool lgMustMalloc = true;~~
85
86	263	DEBUG_ENTRY( "mole_newton_step()" );
87	264
88		if( lgMustMalloc )
89		{
90		/* on very first call must create space */
91		lgMustMalloc = false;
92
93		amat = ((double)MALLOC( (size_t)(mole.num_compactedmole.num_compacted)*sizeof(double) ));
94		b0vec = ((double)MALLOC( (size_t)(mole.num_compacted)sizeof(double) ));
95		b1vec = ((double)MALLOC( (size_t)(mole.num_compacted)sizeof(double) ));
96		b2vec = ((double)MALLOC( (size_t)(mole.num_compacted)sizeof(double) ));
97		calcv = ((double)MALLOC( (size_t)(mole.num_total)sizeof(double) ));
98		b = (double )MALLOC((size_t)mole.num_totalsizeof(double));
99		c = (double *)MALLOC((size_t)mole.num_totalsizeof(double *));
100		c[0] = (double )MALLOC((size_t)mole.num_total
101		mole.num_total*sizeof(double));
102		for(i=1;i<mole.num_total;i++)
103		{
104		c[i] = c[i-1]+mole.num_total;
105		}
106		}
107
108		if( iteration >= dynamics.n_initial_relax+1 && dynamics.lgAdvection /&& radius.depth < dynamics.oldFullDepth/ )
109		{
110		ASSERT(dynamics.Rate > 0.0);
111		lgConserve = false;
112		}
113		else
114		{
115		lgConserve = true;
116		}
117
118		for(i=0;i<mole.num_total;i++)
119		{
120		calcv[i] = mole.list[i]->den;
121		}
122
123		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
124		{
125		if (element_list[nelem]->ipMl[0] != -1)
126		{
127		sum = 0.;
128		for (ion=0; ion<N_MOLE_ION; ion++)
129		{
130		if (element_list[nelem]->ipMl[ion] != -1)
131		sum += calcv[element_list[nelem]->ipMl[ion]];
132		}
133		if (sum > SMALLFLOAT)
134		{
135		for (ion=0; ion<N_MOLE_ION; ion++)
136		{
137		if (element_list[nelem]->ipMl[ion] != -1)
138		fion[nelem][ion] = calcv[element_list[nelem]->ipMl[ion]]/sum;
139		else
140		fion[nelem][ion] = 0.;
141		}
142		}
143		else
144		{
145		lgSet = false;
146		for (ion=0; ion<N_MOLE_ION; ion++)
147		{
148		if (element_list[nelem]->ipMl[ion] != -1 && !lgSet)
149		{
150		fion[nelem][ion] = 1.0;
151		lgSet = true;
152		}
153		else
154		{
155		fion[nelem][ion] = 0.;
156		}
157		}
158		}
159
160		lgSet = false;
161		for (ion=0;ion<N_MOLE_ION;ion++)
162		{
163		if (element_list[nelem]->ipMl[ion] != -1)
164		{
165		if (!lgSet)
166		calcv[element_list[nelem]->ipMl[ion]] = sum;
167		else
168		calcv[element_list[nelem]->ipMl[ion]] = 0.;
169		lgSet = true;
170		}
171		}
172		}
173		}
174
175		for( i=0; i < mole.num_compacted; i++ )
176		{
177		b0vec[i] = b2vec[i] = calcv[groupspecies[i]->index];
	265	MoleMap.setup(&b0vec[0]);
	266
	267	for( i=0; i < n; i++ )
	268	{
	269	b2vec[i] = b0vec[i];
178	270	}
179	271
…	…
187	279	flast = fstep;
188	280
189		/* Generate chemical balance vector (mole.b[]) and Jacobian array
190		(c[][], first iteration) from reaction list */
191
192		mole_eval_dynamic_balance(mole.num_total,b,(loop==0)?c:NULL);
193
194		/* last test - do not include advection if we have overrun the radius scale
195		* of previous iteration */
196
197		/* add positive ions and neutral atoms: ratios are set by ion_solver,
198		* we determine abundance of the group as a whole here */
199
200		if (loop == 0) {
201		// c[i][i] can be positive for auto-catalytic reactions
202		//for(i=0;i<mole.num_calc;i++)
203		//{
204		// ASSERT (c[i][i] <= 0.);
205		//}
206		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
207		{
208		if (element_list[nelem]->ipMl[0] != -1)
209		{
210		for(i=0;i<mole.num_calc;i++)
211		{
212		ASSERT(mole.list[i]->index == i);
213		sum = 0.;
214		for (ion=0;ion<N_MOLE_ION;ion++)
215		{
216		if (element_list[nelem]->ipMl[ion] != -1)
217		{
218		sum += fion[nelem][ion]*c[element_list[nelem]->ipMl[ion]][i];
219		c[element_list[nelem]->ipMl[ion]][i] = 0.;
220		}
221		}
222		c[element_list[nelem]->ipMl[0]][i] = sum;
223		}
224		}
225		}
226		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
227		{
228		if (element_list[nelem]->ipMl[0] != -1)
229		{
230		for(i=0;i<mole.num_calc;i++)
231		{
232		sum = 0.0;
233		for (ion=0;ion<N_MOLE_ION;ion++)
234		{
235		if (element_list[nelem]->ipMl[ion] != -1)
236		{
237		sum += c[i][element_list[nelem]->ipMl[ion]];
238		c[i][element_list[nelem]->ipMl[ion]] = 0.;
239		}
240		}
241		c[i][element_list[nelem]->ipMl[0]] = sum;
242		}
243		}
244		}
245		}
246
247		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
248		{
249		if (element_list[nelem]->ipMl[0] != -1)
250		{
251		sum = 0.0;
252		for (ion=0;ion<N_MOLE_ION;ion++)
253		{
254		if (element_list[nelem]->ipMl[ion] != -1)
255		{
256		sum += b[element_list[nelem]->ipMl[ion]];
257		b[element_list[nelem]->ipMl[ion]] = 0.0;
258		}
259		}
260		b[element_list[nelem]->ipMl[0]] = sum;
261		}
262		}
263
264		/* For Newton-Raphson method, want the change in populations to be zero,
265		* so the conserved component must also be zero */
266		if (lgConserve)
267		{
268		double scale;
269		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
270		{
271		if (element_list[nelem]->ipMl[0] != -1)
272		{
273		if (loop == 0)
274		{
275		scale = -(ABSLIM+fabs(c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]]));
276		for(i=0;i<mole.num_calc;i++)
277		{
278		c[i][element_list[nelem]->ipMl[0]] = mole.list[i]->nElem[nelem]*scale;
279		}
280		}
281		else
282		scale = c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]];
283
284		b[element_list[nelem]->ipMl[0]] = 0.;
285		}
286		}
287		}
	281	funjac(&ervals[0],&amat[0],loop);
	282
	283	if (loop == 0)
	284	{
	285	for( i=0; i < n; i++ )
	286	{
	287	escale[i] = MAT(amat,i,i);
	288	}
	289	}
288	290
289	291	error0 = emax = 0.;
290	292	iworst = -1;
291		for( i=0; i < ~~mole.num_compacted~~; i++ )
	293	for( i=0; i < n; i++ )
292	294	{
293	295	/* Smooth the error mode tailoff */
294		etmp = b[groupspecies[i]->index]/
295		((ERRLIM+fabs(b0vec[i]))*(ERRLIM+fabs(c[groupspecies[i]->index][groupspecies[i]->index])));
	296	etmp = ervals[i]/((ERRLIM+fabs(b0vec[i]))*(ERRLIM+fabs(escale[i])));
296	297	etmp *= etmp;
297	298	if (etmp > emax)
…	…
309	310	// loop,error0,fstep,groupspecies[iworst]->label,groupspecies[iworst]->den,,emax);
310	311
311		~~/------------------------------------------------------------------ /~~
312		~~if(printsol \|\| (trace.lgTrace && trace.lgTr_H2_Mole ))~~
313		{
314		~~/* print the full matrix */~~
315		~~fprintf( ioQQQ, " ");~~
316		~~for( i=0; i < mole.num_calc; i++ )~~
317		{
318		~~fprintf( ioQQQ, " %-4.4s", mole.list[i]->label );~~
319		}
320		~~fprintf( ioQQQ, " \n" );~~
321
322		~~fprintf(ioQQQ," MOLE old abundances\t%.2f",fnzone);~~
323		~~for( i=0; i<mole.num_calc; i++ )~~
324		~~fprintf(ioQQQ,"\t%.2e", mole.list[i]->den );~~
325		~~fprintf(ioQQQ,"\n" );~~
326
327		~~for( i=0; i < mole.num_calc; i++ )~~
328		{
329		~~fprintf( ioQQQ, " MOLE%2ld %-4.4s", i ,mole.list[i]->label);~~
330		~~for( j=0; j < mole.num_calc; j++ )~~
331		{
332		~~fprintf( ioQQQ, "%10.2e", c[j][i] );~~
333		}
334		~~fprintf( ioQQQ, "%10.2e", b[i] );~~
335		~~fprintf( ioQQQ, "\n" );~~
336		}
337		}
338
339	312	if (loop == 0)
340	313	{
341		/* Find descent direction */
342		for( i=0; i < mole.num_compacted; i++ )
343		{
344		b1vec[i] = b[groupspecies[i]->index];
345		// c[i][i] can be +ve for (effectively) autocatalytic reactions,
346		// e.g. H-,H+=>H,H
347		//ASSERT (c[groupspecies[i]->index][groupspecies[i]->index] <= 0.);
348		for( j=0; j < mole.num_compacted; j++ )
349		{
350		MAT(amat,i,j) = c[groupspecies[i]->index][groupspecies[j]->index];
351		}
352		}
353
354
	314	for( i=0; i < n; i++ )
	315	{
	316	b1vec[i] = ervals[i];
	317	}
	318	merror = solve_system(&amat[0],&b1vec[0],n);
355	319	error1 = error0;
356	320	grad = -2*error1;
357		~~merror = solve_system((double *)amat,b1vec,mole.num_compacted);~~
358
359	321	/* b1vec is now descent direction */
360	322
…	…
393	355	if (0)
394	356	{
395		for( i=0; i < ~~mole.num_compacted~~; i++ )
	357	for( i=0; i < n; i++ )
396	358	{
397	359	fprintf(ioQQQ,"%15.7g\t",b1vec[i]);
…	…
403	365	{
404	366	if (1)
405		for (i=0;i<~~mole.num_compacted~~;i++)
	367	for (i=0;i<n;i++)
406	368	{
407	369	if (groupspecies[i] == findspecies("HNC"))
…	…	</