paso/src/SolverFCT.c

/* $Id:$ */

/*******************************************************
 *
 *       Copyright 2007 by University of Queensland
 *
 *                http://esscc.uq.edu.au
 *        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);
           MEMFREE(in->u);
           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;
     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->theta=theta;
     out->dt_max=LARGE_POSITIVE_FLOAT;
     out->valid_matrices=FALSE;
     out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
     Paso_SystemMatrix_allocBuffer(out->transport_matrix);
     out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
     out->iteration_matrix=NULL;

     out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
     out->u=NULL;
     out->u_min=0.;
     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);

         out->u=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);

         if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) || 
                 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix))  ) {
             
             #pragma omp parallel 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 parallel for schedule(static) private(i,iptr,iptr_main)
             for (i = 0; i < n; ++i) {
                    iptr_main=pattern->mainPattern->ptr[0]-1;
                    for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
                          if (pattern->mainPattern->index[iptr]==i) {
                               iptr_main=iptr;
                               break;
                          }
                    }
                    out->main_iptr[i]=iptr_main;
                    if (iptr_main==pattern->mainPattern->ptr[0]-1)
                         Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
             }

      }

  }
  if (Paso_noError()) {
     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_FCTransportProblem_getTotalNumRows(in);
    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) 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
    in->u_min=u_min;
    #pragma omp parallel for schedule(static) private(i)
    for (i=0;i<n;++i) {
        in->u[i]=u[i]-u_min;
    }
}
1	gross	1361	/* $Id:$ */
2
3			/*******************************************************
4			*
5			* Copyright 2007 by University of Queensland
6			*
7			* http://esscc.uq.edu.au
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	gross	1362	/* Paso: FCTransportProblem */
17	gross	1361
18			/**************************************************************/
19
20			/* Author: l.gross@uq.edu.au */
21
22			/**************************************************************/
23
24
25			#include "Paso.h"
26			#include "SolverFCT.h"
27			#include "PasoUtil.h"
28
29
30	gross	1362	/**************************************************************/
31	gross	1361
32	gross	1362	/* free all memory used by */
33			void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {
34			if (in!=NULL) {
35			in->reference_counter--;
36			if (in->reference_counter<=0) {
37	gross	1361
38	gross	1362	Paso_SystemMatrix_free(in->transport_matrix);
39	gross	1407	Paso_SystemMatrix_free(in->mass_matrix);
40			Paso_SystemMatrix_free(in->iteration_matrix);
41	gross	1362	Paso_MPIInfo_free(in->mpi_info);
42	gross	1364	MEMFREE(in->u);
43	gross	1407	MEMFREE(in->main_iptr);
44	gross	1362	MEMFREE(in->lumped_mass_matrix);
45	gross	1407	MEMFREE(in->main_diagonal_low_order_transport_matrix);
46	gross	1362	MEMFREE(in);
47			}
48			}
49			}
50	gross	1361
51	gross	1362	Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
52	gross	1361	if (in!=NULL) {
53	gross	1362	++(in->reference_counter);
54	gross	1361	}
55	gross	1367	return in;
56	gross	1362	}
57
58			Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
59			return in->transport_matrix;
60	gross	1361	}
61	gross	1407	Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
62			return in->mass_matrix;
63	gross	1362	}
64	gross	1361
65	gross	1362	double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
66			return in->lumped_mass_matrix;
67			}
68	gross	1361
69	gross	1362	dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
70			return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
71			}
72	gross	1361
73	gross	1407	Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size
74	gross	1361
75	gross	1362
76			) {
77			Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
78			Paso_FCTransportProblem* out=NULL;
79			dim_t n,i;
80	trankine	1462	index_t iptr,iptr_main;
81	gross	1362
82			if ((theta<0.) \|\| (theta >1.)) {
83			Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
84	gross	1361	return NULL;
85	gross	1362	}
86	gross	1361
87	gross	1362	out=MEMALLOC(1,Paso_FCTransportProblem);
88			if (Paso_checkPtr(out)) return NULL;
89	gross	1361
90	gross	1407
91	gross	1362	out->theta=theta;
92	gross	1407	out->dt_max=LARGE_POSITIVE_FLOAT;
93	gross	1363	out->valid_matrices=FALSE;
94	gross	1362	out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
95			Paso_SystemMatrix_allocBuffer(out->transport_matrix);
96	gross	1407	out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
97			out->iteration_matrix=NULL;
98
99	gross	1362	out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
100	gross	1407	out->u=NULL;
101			out->u_min=0.;
102	gross	1362	out->main_iptr=NULL;
103			out->lumped_mass_matrix=NULL;
104	gross	1407	out->main_diagonal_low_order_transport_matrix=NULL;
105	gross	1361
106	gross	1362	if (Paso_noError()) {
107			n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
108
109	gross	1407	out->u=MEMALLOC(n,double);
110	gross	1362	out->main_iptr=MEMALLOC(n,index_t);
111			out->lumped_mass_matrix=MEMALLOC(n,double);
112	gross	1407	out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
113	gross	1362
114	gross	1407	if ( ! (Paso_checkPtr(out->u) \|\| Paso_checkPtr(out->main_iptr) \|\|
115			Paso_checkPtr(out->lumped_mass_matrix) \|\| Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) ) {
116	gross	1362
117	gross	1367	#pragma omp parallel for schedule(static) private(i)
118			for (i = 0; i < n; ++i) {
119			out->lumped_mass_matrix[i]=0.;
120	gross	1407	out->main_diagonal_low_order_transport_matrix[i]=0.;
121	gross	1367	out->u[i]=0.;
122			}
123	gross	1362	/* identify the main diagonals */
124	ksteube	1558	#pragma omp parallel for schedule(static) private(i,iptr,iptr_main)
125	gross	1362	for (i = 0; i < n; ++i) {
126			iptr_main=pattern->mainPattern->ptr[0]-1;
127			for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
128			if (pattern->mainPattern->index[iptr]==i) {
129			iptr_main=iptr;
130			break;
131			}
132			}
133			out->main_iptr[i]=iptr_main;
134			if (iptr_main==pattern->mainPattern->ptr[0]-1)
135			Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
136			}
137
138			}
139
140	gross	1361	}
141			if (Paso_noError()) {
142			return out;
143			} else {
144	gross	1362	Paso_FCTransportProblem_free(out);
145	gross	1361	return NULL;
146			}
147			}
148	gross	1364
149			void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
150			dim_t i, n;
151	gross	1407	double local_u_min,u_min;
152	gross	1364
153			n=Paso_FCTransportProblem_getTotalNumRows(in);
154	gross	1407	u_min=LARGE_POSITIVE_FLOAT;
155			#pragma omp parallel private(local_u_min)
156			{
157			local_u_min=0.;
158			#pragma omp for schedule(static) private(i)
159	gross	1410	for (i=0;i<n;++i) local_u_min=MIN(local_u_min,u[i]);
160	gross	1407	#pragma omp critical
161			{
162			u_min=MIN(u_min,local_u_min);
163			}
164	gross	1364	}
165	gross	1407	#ifdef PASO_MPI
166			local_u_min=u_min;
167	ksteube	1415	MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
168	gross	1407	#endif
169			in->u_min=u_min;
170	gross	1410	#pragma omp parallel for schedule(static) private(i)
171			for (i=0;i<n;++i) {
172			in->u[i]=u[i]-u_min;
173			}
174	gross	1364	}