paso/src/SolverFCT.c


/*******************************************************
*
* Copyright (c) 2003-2008 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


/**************************************************************/

/* Paso: FCTransportProblem                                          */

/**************************************************************/

/* Author: l.gross@uq.edu.au                                */

/**************************************************************/


#include "Paso.h"
#include "SolverFCT.h"
#include "PasoUtil.h"


/**************************************************************/

/* free all memory used by                                */

void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {
     if (in!=NULL) {
        in->reference_counter--;
        if (in->reference_counter<=0) {
           Paso_SystemMatrix_free(in->transport_matrix);
           Paso_SystemMatrix_free(in->mass_matrix);
           Paso_SystemMatrix_free(in->iteration_matrix);
           Paso_MPIInfo_free(in->mpi_info);
           Paso_Coupler_free(in->u_coupler);
           MEMFREE(in->u);
           MEMFREE(in->constraint_weights);
           MEMFREE(in->main_iptr);
           MEMFREE(in->lumped_mass_matrix);
           MEMFREE(in->main_diagonal_low_order_transport_matrix);
           MEMFREE(in);
        }
    }
}

Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
     if (in!=NULL) {
        ++(in->reference_counter);
     }
     return in;
}    

Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
   return in->transport_matrix;
}
Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
   return in->mass_matrix;
}

double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
    return in->lumped_mass_matrix;
}

dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
    return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
}

Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size) 
{
     Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1;  /* at the moment only block size 1 is supported */
     Paso_FCTransportProblem* out=NULL;
     Paso_SystemMatrixPattern *transport_pattern;
     dim_t n,i;
     index_t iptr,iptr_main;

     if ((theta<0.) || (theta >1.)) {
        Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
        return NULL;
     }

     out=MEMALLOC(1,Paso_FCTransportProblem);
     if (Paso_checkPtr(out)) return NULL;
     out->reference_counter=0;
     out->theta=theta;
     out->dt_max=LARGE_POSITIVE_FLOAT;
     out->constraint_factor=sqrt(LARGE_POSITIVE_FLOAT);
     if (out->theta < 1.) {
            out->dt_factor=MIN(1./(1.-out->theta),DT_FACTOR_MAX);
     } else {
            out->dt_factor=DT_FACTOR_MAX;
     }
     out->valid_matrices=FALSE;
     out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
     out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
     out->iteration_matrix=NULL;
     out->u=NULL;
     out->constraint_weights=NULL;
     out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
     out->main_iptr=NULL;
     out->lumped_mass_matrix=NULL;
     out->main_diagonal_low_order_transport_matrix=NULL;

     if (Paso_noError()) {
         n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
         transport_pattern=out->transport_matrix->pattern;
         out->u=MEMALLOC(n,double);
         out->constraint_weights=MEMALLOC(n,double);
         out->main_iptr=MEMALLOC(n,index_t);
         out->lumped_mass_matrix=MEMALLOC(n,double);
         out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
         out->u_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(out),block_size);

         if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->constraint_weights) || 
                 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Paso_noError()  ) {
             
             #pragma omp parallel 
             {
                 #pragma omp for schedule(static) private(i)
                 for (i = 0; i < n; ++i) {
                    out->lumped_mass_matrix[i]=0.;
                    out->main_diagonal_low_order_transport_matrix[i]=0.;
                    out->u[i]=0.;
                 }
                 /* identify the main diagonals */
                 #pragma omp for schedule(static) private(i,iptr,iptr_main)
                 for (i = 0; i < n; ++i) {
                        iptr_main=transport_pattern->mainPattern->ptr[0]-1;
                        for (iptr=transport_pattern->mainPattern->ptr[i];iptr<transport_pattern->mainPattern->ptr[i+1]; iptr++) {
                              if (transport_pattern->mainPattern->index[iptr]==i) {
                                   iptr_main=iptr;
                                   break;
                              }
                        }
                        out->main_iptr[i]=iptr_main;
                        if (iptr_main==transport_pattern->mainPattern->ptr[0]-1)
                             Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
                 }
     
             }
      }

  }
  if (Paso_noError()) {
     out->reference_counter=1;
     return out;
  } else {
     Paso_FCTransportProblem_free(out);
     return NULL;
  }
} 

void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
    dim_t i, n;
    double local_u_min,u_min;
    n=Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
    u_min=LARGE_POSITIVE_FLOAT;
    #pragma omp parallel private(local_u_min)
    {
         local_u_min=0.;
         #pragma omp for schedule(static) private(i)
         for (i=0;i<n;++i) {
              in->u[i]=u[i];
              local_u_min=MIN(local_u_min,u[i]);
         }
         #pragma omp critical 
         {
            u_min=MIN(u_min,local_u_min);
         }
    }
    #ifdef PASO_MPI
         local_u_min=u_min;
         MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
    #endif
    if (u_min<0) {
       Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_checkinSolution: initial guess must be non-negative.");
    }
}
dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
{
   return in->transport_matrix->row_block_size;
}

Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
{
   return in->transport_matrix->pattern->col_connector;
}

index_t Paso_FCTransportProblem_getTypeId(const index_t solver,const index_t preconditioner, const index_t package,const  bool_t symmetry) 
{
   return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
}

void Paso_FCTransportProblem_setUpConstraint(Paso_FCTransportProblem* fctp,  const double* q, const double factor)
{
   dim_t i, n;
   register double m, rtmp;
   double factor2= fctp->dt_factor * factor;
   n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);

   if ( fctp->valid_matrices ) {
      Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_insertConstraint: you must not insert a constraint is a valid system.");
      return;
   }
   if (factor<=0) {
      Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_insertConstraint: constraint_factor needs to be positive.");
      return;
   }

   
   #pragma omp for schedule(static) private(i,m,rtmp)
   for (i=0;i<n;++i) {
        if (q[i]>0) {
           m=fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]];
           rtmp=factor2 * (m == 0 ? 1 : m);
           fctp->constraint_weights[i]=rtmp;
           fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]]=m+rtmp;
        } else {
           fctp->constraint_weights[i]=0;
        }
   }
   fctp->constraint_factor=factor;
}

void Paso_FCTransportProblem_insertConstraint(Paso_FCTransportProblem* fctp,  const double* r,  double* source)
{
   dim_t i, n;
   n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);

   #pragma omp for schedule(static) private(i)
   for (i=0;i<n;++i) source[i]+=fctp->constraint_weights[i] * r[i];
}
1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2008 by University of Queensland
5	* Earth Systems Science Computational Center (ESSCC)
6	* http://www.uq.edu.au/esscc
7	*
8	* Primary Business: Queensland, Australia
9	* Licensed under the Open Software License version 3.0
10	* http://www.opensource.org/licenses/osl-3.0.php
11	*
12	*******************************************************/
13
14
15	/**************************************************************/
16
17	/* Paso: FCTransportProblem */
18
19	/**************************************************************/
20
21	/* Author: l.gross@uq.edu.au */
22
23	/**************************************************************/
24
25
26	#include "Paso.h"
27	#include "SolverFCT.h"
28	#include "PasoUtil.h"
29
30
31	/**************************************************************/
32
33	/* free all memory used by */
34
35	void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {
36	if (in!=NULL) {
37	in->reference_counter--;
38	if (in->reference_counter<=0) {
39	Paso_SystemMatrix_free(in->transport_matrix);
40	Paso_SystemMatrix_free(in->mass_matrix);
41	Paso_SystemMatrix_free(in->iteration_matrix);
42	Paso_MPIInfo_free(in->mpi_info);
43	Paso_Coupler_free(in->u_coupler);
44	MEMFREE(in->u);
45	MEMFREE(in->constraint_weights);
46	MEMFREE(in->main_iptr);
47	MEMFREE(in->lumped_mass_matrix);
48	MEMFREE(in->main_diagonal_low_order_transport_matrix);
49	MEMFREE(in);
50	}
51	}
52	}
53
54	Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
55	if (in!=NULL) {
56	++(in->reference_counter);
57	}
58	return in;
59	}
60
61	Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
62	return in->transport_matrix;
63	}
64	Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
65	return in->mass_matrix;
66	}
67
68	double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
69	return in->lumped_mass_matrix;
70	}
71
72	dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
73	return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
74	}
75
76	Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size)
77	{
78	Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
79	Paso_FCTransportProblem* out=NULL;
80	Paso_SystemMatrixPattern *transport_pattern;
81	dim_t n,i;
82	index_t iptr,iptr_main;
83
84	if ((theta<0.) \|\| (theta >1.)) {
85	Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
86	return NULL;
87	}
88
89	out=MEMALLOC(1,Paso_FCTransportProblem);
90	if (Paso_checkPtr(out)) return NULL;
91	out->reference_counter=0;
92	out->theta=theta;
93	out->dt_max=LARGE_POSITIVE_FLOAT;
94	out->constraint_factor=sqrt(LARGE_POSITIVE_FLOAT);
95	if (out->theta < 1.) {
96	out->dt_factor=MIN(1./(1.-out->theta),DT_FACTOR_MAX);
97	} else {
98	out->dt_factor=DT_FACTOR_MAX;
99	}
100	out->valid_matrices=FALSE;
101	out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
102	out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
103	out->iteration_matrix=NULL;
104	out->u=NULL;
105	out->constraint_weights=NULL;
106	out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
107	out->main_iptr=NULL;
108	out->lumped_mass_matrix=NULL;
109	out->main_diagonal_low_order_transport_matrix=NULL;
110
111	if (Paso_noError()) {
112	n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
113	transport_pattern=out->transport_matrix->pattern;
114	out->u=MEMALLOC(n,double);
115	out->constraint_weights=MEMALLOC(n,double);
116	out->main_iptr=MEMALLOC(n,index_t);
117	out->lumped_mass_matrix=MEMALLOC(n,double);
118	out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
119	out->u_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(out),block_size);
120
121	if ( ! (Paso_checkPtr(out->u) \|\| Paso_checkPtr(out->main_iptr) \|\| Paso_checkPtr(out->constraint_weights) \|\|
122	Paso_checkPtr(out->lumped_mass_matrix) \|\| Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Paso_noError() ) {
123
124	#pragma omp parallel
125	{
126	#pragma omp for schedule(static) private(i)
127	for (i = 0; i < n; ++i) {
128	out->lumped_mass_matrix[i]=0.;
129	out->main_diagonal_low_order_transport_matrix[i]=0.;
130	out->u[i]=0.;
131	}
132	/* identify the main diagonals */
133	#pragma omp for schedule(static) private(i,iptr,iptr_main)
134	for (i = 0; i < n; ++i) {
135	iptr_main=transport_pattern->mainPattern->ptr[0]-1;
136	for (iptr=transport_pattern->mainPattern->ptr[i];iptr<transport_pattern->mainPattern->ptr[i+1]; iptr++) {
137	if (transport_pattern->mainPattern->index[iptr]==i) {
138	iptr_main=iptr;
139	break;
140	}
141	}
142	out->main_iptr[i]=iptr_main;
143	if (iptr_main==transport_pattern->mainPattern->ptr[0]-1)
144	Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
145	}
146
147	}
148	}
149
150	}
151	if (Paso_noError()) {
152	out->reference_counter=1;
153	return out;
154	} else {
155	Paso_FCTransportProblem_free(out);
156	return NULL;
157	}
158	}
159
160	void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
161	dim_t i, n;
162	double local_u_min,u_min;
163	n=Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
164	u_min=LARGE_POSITIVE_FLOAT;
165	#pragma omp parallel private(local_u_min)
166	{
167	local_u_min=0.;
168	#pragma omp for schedule(static) private(i)
169	for (i=0;i<n;++i) {
170	in->u[i]=u[i];
171	local_u_min=MIN(local_u_min,u[i]);
172	}
173	#pragma omp critical
174	{
175	u_min=MIN(u_min,local_u_min);
176	}
177	}
178	#ifdef PASO_MPI
179	local_u_min=u_min;
180	MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
181	#endif
182	if (u_min<0) {
183	Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_checkinSolution: initial guess must be non-negative.");
184	}
185	}
186	dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
187	{
188	return in->transport_matrix->row_block_size;
189	}
190
191	Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
192	{
193	return in->transport_matrix->pattern->col_connector;
194	}
195
196	index_t Paso_FCTransportProblem_getTypeId(const index_t solver,const index_t preconditioner, const index_t package,const bool_t symmetry)
197	{
198	return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
199	}
200
201	void Paso_FCTransportProblem_setUpConstraint(Paso_FCTransportProblem* fctp, const double* q, const double factor)
202	{
203	dim_t i, n;
204	register double m, rtmp;
205	double factor2= fctp->dt_factor * factor;
206	n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
207
208	if ( fctp->valid_matrices ) {
209	Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_insertConstraint: you must not insert a constraint is a valid system.");
210	return;
211	}
212	if (factor<=0) {
213	Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_insertConstraint: constraint_factor needs to be positive.");
214	return;
215	}
216
217
218	#pragma omp for schedule(static) private(i,m,rtmp)
219	for (i=0;i<n;++i) {
220	if (q[i]>0) {
221	m=fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]];
222	rtmp=factor2 * (m == 0 ? 1 : m);
223	fctp->constraint_weights[i]=rtmp;
224	fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]]=m+rtmp;
225	} else {
226	fctp->constraint_weights[i]=0;
227	}
228	}
229	fctp->constraint_factor=factor;
230	}
231
232	void Paso_FCTransportProblem_insertConstraint(Paso_FCTransportProblem* fctp, const double* r, double* source)
233	{
234	dim_t i, n;
235	n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
236
237	#pragma omp for schedule(static) private(i)
238	for (i=0;i<n;++i) source[i]+=fctp->constraint_weights[i] * r[i];
239	}