Context Navigation

← Previous Change
Next Change →

mole_solve.cpp

Timestamp:

05/13/08 17:01:00 (8 months ago)

Author:

rjrw

Message:

Move further towards Householder projection method for of conservation
constraints

Some other tidyings (move from explicit pointers to multi_arr)

Files:

: 1 modified

branches/newmole/source/mole_solve.cpp (modified) (7 diffs)

Legend:

: Unmodified
: Added
: Removed

branches/newmole/source/mole_solve.cpp

r2054	r2055
308	308
309	309
310		STATIC void mole_eval_dynamic_balance(long int num_total, double b, double *c);
	310	STATIC void mole_eval_dynamic_balance(long int num_total, double *b, bool lgJac, multi_arr<double,2> &c);
311	311	enum {PRINTSOL = false};
312	312
…	…
339	339	j;
340	340	double sum;
341		static double
342		**c;
	341	multi_arr<double,2> c(mole.num_total,mole.num_total);
343	342	int printsol = PRINTSOL;
344	343	valarray<double> b(mole.num_total);
345		~~static bool lgMustMalloc = true;~~
346	344	bool lgConserve;
347	345
…	…
360	358	}
361	359
362		~~if( lgMustMalloc )~~
363		{
364		~~/* on very first call must create space */~~
365		~~lgMustMalloc = false;~~
366		~~c = (double *)MALLOC((size_t)mole.num_totalsizeof(double *));~~
367		c[0] = (double )MALLOC((size_t)mole.num_total
368		~~mole.num_total*sizeof(double));~~
369		~~for(i=1;i<mole.num_total;i++)~~
370		{
371		~~c[i] = c[i-1]+mole.num_total;~~
372		}
373		}
374
375	360	/* Generate chemical balance vector (mole.b[]) and Jacobian array
376	361	(c[][], first iteration) from reaction list */
377	362
378		mole_eval_dynamic_balance(mole.num_total,&b[0],lgJac?c:NULL);
379
380		/* add positive ions and neutral atoms: ratios are set by ion_solver,
381		* we determine abundance of the group as a whole here */
382
383		if (lgJac) {
384		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
385		{
386		if (element_list[nelem]->ipMl[0] != -1)
387		{
388		for(i=0;i<mole.num_calc;i++)
389		{
390		ASSERT(mole.list[i]->index == i);
391		sum = 0.;
392		for (ion=0;ion<N_MOLE_ION;ion++)
393		{
394		if (element_list[nelem]->ipMl[ion] != -1)
395		{
396		sum += MoleMap.fion[nelem][ion]*c[element_list[nelem]->ipMl[ion]][i];
397		c[element_list[nelem]->ipMl[ion]][i] = 0.;
398		}
399		}
400		c[element_list[nelem]->ipMl[0]][i] = sum;
401		}
402		}
403		}
404		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
405		{
406		if (element_list[nelem]->ipMl[0] != -1)
407		{
408		for(i=0;i<mole.num_calc;i++)
409		{
410		sum = 0.0;
411		for (ion=0;ion<N_MOLE_ION;ion++)
412		{
413		if (element_list[nelem]->ipMl[ion] != -1)
414		{
415		sum += c[i][element_list[nelem]->ipMl[ion]];
416		c[i][element_list[nelem]->ipMl[ion]] = 0.;
417		}
418		}
419		c[i][element_list[nelem]->ipMl[0]] = sum;
420		}
421		}
422		}
423		}
424
425		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
426		{
427		if (element_list[nelem]->ipMl[0] != -1)
428		{
429		sum = 0.0;
430		for (ion=0;ion<N_MOLE_ION;ion++)
431		{
432		if (element_list[nelem]->ipMl[ion] != -1)
433		{
434		sum += b[element_list[nelem]->ipMl[ion]];
435		b[element_list[nelem]->ipMl[ion]] = 0.0;
436		}
437		}
438		b[element_list[nelem]->ipMl[0]] = sum;
439		}
440		}
441
442		/* For Newton-Raphson method, want the change in populations to be zero,
443		* so the conserved component must also be zero */
444		if (lgConserve)
445		{
446		// Development code towards using Householder reflections to apply
447		// conservation constraints
448		if (0)
449		{
450		multi_arr<double,2> u(mole.num_calc,LIMELM);
451		valarray<double> betah(LIMELM);
452		for(i=0;i<mole.num_calc;i++)
453		{
454		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
455		{
456		u[i][nelem] = mole.list[i]->nElem[nelem];
457		}
458		}
459		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
460		{
461		double beta;
462		valarray<double> x(mole.num_calc-nelem),
463		y(mole.num_calc-nelem),
464		v(mole.num_calc-nelem),
465		p(mole.num_calc-nelem);
466		for(i=0;i<mole.num_calc-nelem;i++)
467		{
468		x[i] = u[i+nelem][nelem];
469		}
470		householder(x,v,beta);
471		betah[nelem] = beta;
472		for(i=0;i<mole.num_calc-nelem;i++)
473		{
474		u[i+nelem][nelem] = v[i];
475		}
476		for(long nother=nelem+1;nother<LIMELM;nother++)
477		{
478		for (i=0;i<mole.num_calc-nelem;i++)
479		{
480		y[i] = u[i+nelem][nother];
481		}
482		p = project(y,v,beta);
483		for (i=0;i<mole.num_calc-nelem;i++)
484		{
485		u[i+nelem][nother] = p[i];
486		}
487		}
488		}
489
490		valarray<double> b2(mole.num_calc);
491		for(i=0;i<mole.num_calc;i++)
492		{
493		b2[i] = b[i];
494		}
495		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
496		{
497		valarray<double> v(mole.num_calc-nelem),
498		b3(mole.num_calc-nelem),
499		p(mole.num_calc-nelem);
500		for(i=0;i<mole.num_calc-nelem;i++)
501		{
502		v[i] = u[i+nelem][nelem];
503		b3[i] = b2[i+nelem];
504		}
505		p = project(b3,v,betah[nelem]);
506		for (i=0;i<mole.num_calc-nelem;i++)
507		{
508		b2[i+nelem] = p[i];
509		}
510		}
511		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
512		{
513		fprintf(ioQQQ,"Elem %ld proj %g %g\n",nelem,b2[nelem],dense.gas_phase[nelem]);
514		}
515		for (;nelem<mole.num_calc;nelem++)
516		{
517		fprintf(ioQQQ,"Remain %ld proj %g\n",nelem,b2[nelem]);
518		}
519
520		for(long j=0;j<mole.num_calc;++j)
521		{
522		for(i=0;i<mole.num_calc;++i)
523		{
524		b2[i] = c[j][i];
525		}
526		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
527		{
528		valarray<double> v(mole.num_calc-nelem),
529		b3(mole.num_calc-nelem),
530		p(mole.num_calc-nelem);
531		for(i=0;i<mole.num_calc-nelem;i++)
532		{
533		v[i] = u[i+nelem][nelem];
534		b3[i] = b2[i+nelem];
535		}
536		p = project(b3,v,betah[nelem]);
537		for (i=0;i<mole.num_calc-nelem;i++)
538		{
539		b2[i+nelem] = p[i];
540		}
541		}
542		fprintf(ioQQQ,"Eqn... %ld\n",j);
543		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
544		{
545		fprintf(ioQQQ,"Eqqq %ld proj %g\n",nelem,b2[nelem]);
546		}
547		for (;nelem<mole.num_calc;nelem++)
548		{
549		fprintf(ioQQQ,"Remn %ld proj %g\n",nelem,b2[nelem]);
550		}
551		}
552		}
553		if (0)
554		{
555		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
556		{
557		double sum1 = 0.;
558		for(i=0;i<mole.num_calc;i++)
559		{
560		sum1 += b[i]*mole.list[i]->nElem[nelem];
561		}
562		fprintf(ioQQQ,"Elem %ld cons %g %g\n",nelem,sum1,dense.gas_phase[nelem]);
563		if (lgJac)
564		{
565		for(i=0;i<mole.num_calc;i++)
566		{
567		sum1 = 0.;
568		for (long j=0;j<mole.num_calc;j++)
569		{
570		sum1 += c[i][j]*mole.list[j]->nElem[nelem];
571		}
572		fprintf(ioQQQ,"Eqn %ld error %g\n",i,sum1);
573		}
574		}
575		}
576		}
577
578		double scale;
579		for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
580		{
581		if (element_list[nelem]->ipMl[0] != -1)
582		{
583		if (lgJac)
584		{
585		scale = -(ABSLIM+fabs(c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]]));
586		for(i=0;i<mole.num_calc;i++)
587		{
588		c[i][element_list[nelem]->ipMl[0]] = mole.list[i]->nElem[nelem]*scale;
589		}
590		}
591		else
592		{
593		scale = c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]];
594		}
595
596		b[element_list[nelem]->ipMl[0]] = 0.;
597		}
598		}
599		}
600
	363	mole_eval_dynamic_balance(mole.num_total,&b[0],lgJac,c);
	364
601	365	/------------------------------------------------------------------ /
602	366	if(printsol \|\| (trace.lgTrace && trace.lgTr_H2_Mole ))
…	…
627	391	}
628	392
629		for( i=0; i < mole.num_compacted; i++ )
630		{
631		ervals[i] = b[groupspecies[i]->index];
632		}
633
634		if (lgJac)
635		{
	393	/* add positive ions and neutral atoms: ratios are set by ion_solver,
	394	* we determine abundance of the group as a whole here */
	395
	396	if (lgJac) {
	397	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	398	{
	399	if (element_list[nelem]->ipMl[0] != -1)
	400	{
	401	for(i=0;i<mole.num_calc;i++)
	402	{
	403	ASSERT(mole.list[i]->index == i);
	404	sum = 0.;
	405	for (ion=0;ion<N_MOLE_ION;ion++)
	406	{
	407	if (element_list[nelem]->ipMl[ion] != -1)
	408	{
	409	sum += MoleMap.fion[nelem][ion]*c[element_list[nelem]->ipMl[ion]][i];
	410	c[element_list[nelem]->ipMl[ion]][i] = 0.;
	411	}
	412	}
	413	c[element_list[nelem]->ipMl[0]][i] = sum;
	414	}
	415	}
	416	}
	417	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	418	{
	419	if (element_list[nelem]->ipMl[0] != -1)
	420	{
	421	for(i=0;i<mole.num_calc;i++)
	422	{
	423	sum = 0.0;
	424	for (ion=0;ion<N_MOLE_ION;ion++)
	425	{
	426	if (element_list[nelem]->ipMl[ion] != -1)
	427	{
	428	sum += c[i][element_list[nelem]->ipMl[ion]];
	429	c[i][element_list[nelem]->ipMl[ion]] = 0.;
	430	}
	431	}
	432	c[i][element_list[nelem]->ipMl[0]] = sum;
	433	}
	434	}
	435	}
	436	}
	437
	438	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	439	{
	440	if (element_list[nelem]->ipMl[0] != -1)
	441	{
	442	sum = 0.0;
	443	for (ion=0;ion<N_MOLE_ION;ion++)
	444	{
	445	if (element_list[nelem]->ipMl[ion] != -1)
	446	{
	447	sum += b[element_list[nelem]->ipMl[ion]];
	448	b[element_list[nelem]->ipMl[ion]] = 0.0;
	449	}
	450	}
	451	b[element_list[nelem]->ipMl[0]] = sum;
	452	}
	453	}
	454
	455	/* Species are now grouped -- only mole.num_compacted elements now active */
	456
	457	/* For Newton-Raphson method, want the change in populations to be zero,
	458	* so the conserved component must also be zero */
	459	if (lgConserve)
	460	{
	461	if (0)
	462	{
	463	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	464	{
	465	double sum1 = 0.;
	466	for(i=0;i<mole.num_calc;i++)
	467	{
	468	sum1 += b[i]*mole.list[i]->nElem[nelem];
	469	}
	470	fprintf(ioQQQ,"Elem %ld cons %g %g\n",nelem,sum1,dense.gas_phase[nelem]);
	471	if (lgJac)
	472	{
	473	for(i=0;i<mole.num_calc;i++)
	474	{
	475	sum1 = 0.;
	476	for (long j=0;j<mole.num_calc;j++)
	477	{
	478	sum1 += c[i][j]*mole.list[j]->nElem[nelem];
	479	}
	480	fprintf(ioQQQ,"Eqn %ld error %g\n",i,sum1);
	481	}
	482	}
	483	}
	484	}
	485	}
	486
	487	// Switch to (0) for development code towards using Householder
	488	// reflections to apply conservation constraints
	489	if (1)
	490	{
	491	if (lgConserve)
	492	{
	493	// Original scheme -- replace atom rows with conservation constraints
	494	double scale;
	495	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	496	{
	497	if (element_list[nelem]->ipMl[0] != -1)
	498	{
	499	if (lgJac)
	500	{
	501	scale = -(ABSLIM+fabs(c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]]));
	502	for(i=0;i<mole.num_calc;i++)
	503	{
	504	c[i][element_list[nelem]->ipMl[0]] = mole.list[i]->nElem[nelem]*scale;
	505	}
	506	}
	507	else
	508	{
	509	scale = c[element_list[nelem]->ipMl[0]][element_list[nelem]->ipMl[0]];
	510	}
	511
	512	b[element_list[nelem]->ipMl[0]] = 0.;
	513	}
	514	}
	515	}
	516
636	517	for( i=0; i < mole.num_compacted; i++ )
637	518	{
638		for( j=0; j < mole.num_compacted; j++ )
639		{
640		MAT(amat,i,j) = c[groupspecies[i]->index][groupspecies[j]->index];
641		}
	519	ervals[i] = b[groupspecies[i]->index];
	520	}
	521
	522	if (lgJac)
	523	{
	524	for( i=0; i < mole.num_compacted; i++ )
	525	{
	526	for( j=0; j < mole.num_compacted; j++ )
	527	{
	528	MAT(amat,i,j) = c[groupspecies[i]->index][groupspecies[j]->index];
	529	}
	530	}
	531	}
	532	}
	533	else
	534	{
	535	if (lgConserve)
	536	{
	537	// Prepare Householder vectors for conservation constraints
	538	// -- should be calculated once outside solution loop
	539	multi_arr<double,2> u(mole.num_compacted,LIMELM);
	540	valarray<double> betah(LIMELM);
	541	for(i=0;i<mole.num_compacted;i++)
	542	{
	543	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	544	{
	545	u[i][nelem] = mole.list[groupspecies[i]->index]->nElem[nelem];
	546	}
	547	}
	548	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	549	{
	550	double beta;
	551	valarray<double> x(mole.num_compacted-nelem),
	552	y(mole.num_compacted-nelem),
	553	v(mole.num_compacted-nelem),
	554	p(mole.num_compacted-nelem);
	555	for(i=0;i<mole.num_compacted-nelem;i++)
	556	{
	557	x[i] = u[i+nelem][nelem];
	558	}
	559	householder(x,v,beta);
	560	betah[nelem] = beta;
	561	for(i=0;i<mole.num_compacted-nelem;i++)
	562	{
	563	u[i+nelem][nelem] = v[i];
	564	}
	565	for(long nother=nelem+1;nother<LIMELM;nother++)
	566	{
	567	for (i=0;i<mole.num_compacted-nelem;i++)
	568	{
	569	y[i] = u[i+nelem][nother];
	570	}
	571	p = project(y,v,beta);
	572	for (i=0;i<mole.num_compacted-nelem;i++)
	573	{
	574	u[i+nelem][nother] = p[i];
	575	}
	576	}
	577	}
	578
	579	// Project conservation constraints out of error vector
	580	valarray<double> b2(mole.num_compacted);
	581	for(i=0;i<mole.num_compacted;i++)
	582	{
	583	b2[i] = b[groupspecies[i]->index];
	584	}
	585	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	586	{
	587	valarray<double> v(mole.num_compacted-nelem),
	588	b3(mole.num_compacted-nelem),
	589	p(mole.num_compacted-nelem);
	590	for(i=0;i<mole.num_compacted-nelem;i++)
	591	{
	592	v[i] = u[i+nelem][nelem];
	593	b3[i] = b2[i+nelem];
	594	}
	595	p = project(b3,v,betah[nelem]);
	596	for (i=0;i<mole.num_compacted-nelem;i++)
	597	{
	598	b2[i+nelem] = p[i];
	599	ervals[i] = p[i];
	600	}
	601	}
	602
	603	// Project conservation constraints from Jacobian
	604	valarray<double> bp(mole.num_compacted);
	605	multi_arr<double,2> c2(mole.num_compacted,mole.num_compacted);
	606	for(long j=0;j<mole.num_compacted;++j)
	607	{
	608	for(i=0;i<mole.num_compacted;++i)
	609	{
	610	bp[i] = c[groupspecies[j]->index][groupspecies[i]->index];
	611	}
	612	for (nelem=ipHYDROGEN;nelem<LIMELM;nelem++)
	613	{
	614	valarray<double> v(mole.num_compacted-nelem),
	615	b3(mole.num_compacted-nelem),
	616	p(mole.num_compacted-nelem);
	617	for(i=0;i<mole.num_compacted-nelem;i++)
	618	{
	619	v[i] = u[i+nelem][nelem];
	620	b3[i] = bp[i+nelem];
	621	}
	622	p = project(b3,v,betah[nelem]);
	623	for (i=0;i<mole.num_compacted-nelem;i++)
	624	{
	625	bp[i+nelem] = p[i];
	626	}
	627	}
	628	for (i=0;i<nelem;i++)
	629