/[escript]/trunk/paso/src/ReactiveSolver.c
ViewVC logotype

Contents of /trunk/paso/src/ReactiveSolver.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3005 - (show annotations)
Thu Apr 22 05:59:31 2010 UTC (9 years, 5 months ago) by gross
File MIME type: text/plain
File size: 3858 byte(s)
early call of setPreconditioner in FCT solver removed.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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: reactive solver (D is a diagonal matrix)
18 *
19 * - Mv_t=Dv+q v(0)=u
20 *
21 * to return v(dt)
22 *
23 */
24
25 /**************************************************************/
26
27 /* Author: l.gross@uq.edu.au */
28
29 /**************************************************************/
30
31
32 #include "ReactiveSolver.h"
33 #include "PasoUtil.h"
34 #include "Solver.h"
35
36 err_t Paso_ReactiveSolver_solve(Paso_ReactiveSolver* support, Paso_TransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options, Paso_Performance *pp)
37 {
38 return SOLVER_NO_ERROR;
39 }
40
41 Paso_ReactiveSolver* Paso_ReactiveSolver_alloc(Paso_TransportProblem* fctp)
42 {
43 Paso_ReactiveSolver* out=NULL;
44 out=MEMALLOC(1,Paso_ReactiveSolver);
45 if (Paso_checkPtr(out)) return NULL;
46 return out;
47 }
48
49 void Paso_ReactiveSolver_free(Paso_ReactiveSolver* in)
50 {
51 if (in!=NULL) {
52 MEMFREE(in);
53 }
54 }
55 double Paso_ReactiveSolver_getSafeTimeStepSize(Paso_TransportProblem* fctp)
56 {
57 double dt_max=LARGE_POSITIVE_FLOAT;
58
59
60
61 double beta_0=0., beta_0_loc, betap, beta=EPSILON;
62 dim_t i, n;
63 int fail=0;
64 double dt_max_loc;
65 register double m_ii,m_i, d_ii;
66 #ifdef PASO_MPI
67 double rtmp[2], rtmp_loc[2];
68 #endif
69 n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
70
71 #pragma omp parallel private(beta_0_loc)
72 {
73 beta_0_loc=0.;
74 #pragma omp for schedule(static) private(i,m_ii,m_i, d_ii)
75 for (i=0;i<n;++i) {
76 m_i=fctp->lumped_mass_matrix[i];
77 m_ii=fctp->main_diagonal_mass_matrix[i];
78 /* printf(" %d : %e %e -> %e\n", i, m_i, m_ii, m_i/m_ii); */
79 if ( m_ii>0 ) {
80 beta_0_loc=MAX(m_i/m_ii, beta_0_loc);
81 } else {
82 fail=1;
83 }
84 }
85 #pragma omp critical
86 {
87 beta_0=MAX(beta_0,beta_0_loc);
88 }
89 }
90 #ifdef PASO_MPI
91 rtmp_loc[0]=beta_0;
92 rtmp_loc[1]=1.(*(double) fail );
93 MPI_Allreduce(rtmp_loc, rtmp, 2, MPI_DOUBLE, MPI_MAX, fctp->mpi_info->comm);
94 beta_0=rtmp_loc[0];
95 fail= (rtmp_loc[0] > 0) ? 1 : 0;
96 #endif
97 if (fail>0) {
98 Paso_setError(TYPE_ERROR,"Paso_ReactiveSolver_getSafeTimeStepSize: non-positive main diagonal entry in mass matrix.");
99 }
100
101 if ( Paso_noError()) {
102 if ( ( beta_0 < 2. ) && ( beta_0 >= 1. ) ) {
103 beta=MAX(2*(beta_0-1.)/beta_0, EPSILON);
104 } else {
105 Paso_setError(TYPE_ERROR,"Paso_ReactiveSolver_getSafeTimeStepSize: non-positive main diagonal entry in mass matrix.");
106 }
107 }
108 printf("beta = %e from beta_0 = %e\n",beta, beta_0-1);
109 if ( Paso_noError()) {
110
111 #pragma omp parallel private(dt_max_loc, betap)
112 {
113 betap=beta+1.;
114 dt_max_loc=LARGE_POSITIVE_FLOAT;
115 #pragma omp for schedule(static) private(i, m_ii, m_i, d_ii)
116 for (i=0;i<n;++i) {
117 d_ii=fctp->reactive_matrix[i];
118 m_i=fctp->lumped_mass_matrix[i];
119 m_ii=fctp->main_diagonal_mass_matrix[i];
120
121 if (d_ii > 0.) dt_max_loc=MIN(dt_max_loc, (betap* m_ii - m_i)/d_ii );
122 }
123 #pragma omp critical
124 {
125 dt_max=MIN(dt_max,dt_max_loc);
126 }
127 }
128 dt_max=1./(beta* Paso_Transport_getTheta(fctp));
129 #ifdef PASO_MPI
130 dt_max_loc=dt_max;
131 MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
132 #endif
133 }
134 return dt_max;
135 }

  ViewVC Help
Powered by ViewVC 1.1.26