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);
           Paso_Coupler_free(in->u_coupler);
           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);
     out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
     out->iteration_matrix=NULL;
     out->u=NULL;
     out->u_coupler=Paso_Coupler_alloc(pattern->col_connector,block_size);
     out->u_min=0.;
     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);

         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 
             {
                 #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=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*0;
    #pragma omp parallel for schedule(static) private(i)
    for (i=0;i<n;++i) {
        in->u[i]=u[i]-(in->u_min);
    }
    Paso_Coupler_startCollect(in->u_coupler,in->u);
    Paso_Coupler_finishCollect(in->u_coupler);
}
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;
}
1	/* $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	/* Paso: FCTransportProblem */
17
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	/**************************************************************/
31
32	/* 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	Paso_SystemMatrix_free(in->transport_matrix);
38	Paso_SystemMatrix_free(in->mass_matrix);
39	Paso_SystemMatrix_free(in->iteration_matrix);
40	Paso_MPIInfo_free(in->mpi_info);
41	MEMFREE(in->u);
42	Paso_Coupler_free(in->u_coupler);
43	MEMFREE(in->main_iptr);
44	MEMFREE(in->lumped_mass_matrix);
45	MEMFREE(in->main_diagonal_low_order_transport_matrix);
46	MEMFREE(in);
47	}
48	}
49	}
50
51	Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
52	if (in!=NULL) {
53	++(in->reference_counter);
54	}
55	return in;
56	}
57
58	Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
59	return in->transport_matrix;
60	}
61	Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
62	return in->mass_matrix;
63	}
64
65	double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
66	return in->lumped_mass_matrix;
67	}
68
69	dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
70	return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
71	}
72
73	Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size
74
75
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	index_t iptr,iptr_main;
81
82	if ((theta<0.) \|\| (theta >1.)) {
83	Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
84	return NULL;
85	}
86
87	out=MEMALLOC(1,Paso_FCTransportProblem);
88	if (Paso_checkPtr(out)) return NULL;
89
90
91	out->theta=theta;
92	out->dt_max=LARGE_POSITIVE_FLOAT;
93	out->valid_matrices=FALSE;
94	out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
95	out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
96	out->iteration_matrix=NULL;
97	out->u=NULL;
98	out->u_coupler=Paso_Coupler_alloc(pattern->col_connector,block_size);
99	out->u_min=0.;
100	out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
101	out->main_iptr=NULL;
102	out->lumped_mass_matrix=NULL;
103	out->main_diagonal_low_order_transport_matrix=NULL;
104
105	if (Paso_noError()) {
106	n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
107
108	out->u=MEMALLOC(n,double);
109	out->main_iptr=MEMALLOC(n,index_t);
110	out->lumped_mass_matrix=MEMALLOC(n,double);
111	out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
112
113	if ( ! (Paso_checkPtr(out->u) \|\| Paso_checkPtr(out->main_iptr) \|\|
114	Paso_checkPtr(out->lumped_mass_matrix) \|\| Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) ) {
115
116	#pragma omp parallel
117	{
118	#pragma omp for schedule(static) private(i)
119	for (i = 0; i < n; ++i) {
120	out->lumped_mass_matrix[i]=0.;
121	out->main_diagonal_low_order_transport_matrix[i]=0.;
122	out->u[i]=0.;
123	}
124	/* identify the main diagonals */
125	#pragma omp for schedule(static) private(i,iptr,iptr_main)
126	for (i = 0; i < n; ++i) {
127	iptr_main=pattern->mainPattern->ptr[0]-1;
128	for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
129	if (pattern->mainPattern->index[iptr]==i) {
130	iptr_main=iptr;
131	break;
132	}
133	}
134	out->main_iptr[i]=iptr_main;
135	if (iptr_main==pattern->mainPattern->ptr[0]-1)
136	Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
137	}
138
139	}
140	}
141
142	}
143	if (Paso_noError()) {
144	return out;
145	} else {
146	Paso_FCTransportProblem_free(out);
147	return NULL;
148	}
149	}
150
151	void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
152	dim_t i, n;
153	double local_u_min,u_min;
154	n=Paso_FCTransportProblem_getTotalNumRows(in);
155
156	u_min=LARGE_POSITIVE_FLOAT;
157	#pragma omp parallel private(local_u_min)
158	{
159	local_u_min=0.;
160	#pragma omp for schedule(static) private(i)
161	for (i=0;i<n;++i) local_u_min=MIN(local_u_min,u[i]);
162	#pragma omp critical
163	{
164	u_min=MIN(u_min,local_u_min);
165	}
166	}
167	#ifdef PASO_MPI
168	local_u_min=u_min;
169	MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
170	#endif
171	in->u_min=u_min*0;
172	#pragma omp parallel for schedule(static) private(i)
173	for (i=0;i<n;++i) {
174	in->u[i]=u[i]-(in->u_min);
175	}
176	Paso_Coupler_startCollect(in->u_coupler,in->u);
177	Paso_Coupler_finishCollect(in->u_coupler);
178	}
179	dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
180	{
181	return in->transport_matrix->row_block_size;
182	}
183
184	Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
185	{
186	return in->transport_matrix->pattern->col_connector;
187	}