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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3005 - (show annotations)
Thu Apr 22 05:59:31 2010 UTC (9 years, 6 months ago) by gross
File MIME type: text/plain
File size: 6118 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: Transport solver
18 *
19 * solves Mu_t=Ku+q
20 *
21 * using an operator splitting approach (K=L+D where L is row sum zero and D is a diagonal matrix):
22 *
23 * - Mu_t=Du+q u(0)=u (reactive part)
24 * - Mv_t=Lv v(0)=u(dt/2) (transport part uses flux correction (FCT))
25 * - Mu_t=Du+q u(dt/2)=v(dt/2) (reactive part)
26 *
27 * to return u(dt)
28 *
29 */
30 /**************************************************************/
31
32 /* Author: l.gross@uq.edu.au */
33
34 /**************************************************************/
35
36 #include "Transport.h"
37 #include "FCTSolver.h"
38 #include "ReactiveSolver.h"
39 #include "Solver.h"
40 #include "PasoUtil.h"
41
42
43 void Paso_TransportProblem_solve(Paso_TransportProblem* fctp, double* u, double dt, double* u0, double* q, Paso_Options* options) {
44 Paso_Performance pp;
45 Paso_Function *fctfunction=NULL;
46 Paso_ReactiveSolver *rsolver=NULL;
47 const dim_t FAILURES_MAX=100;
48 dim_t i_substeps, Failed=0;
49 double *u_save=NULL;
50 double dt_max=LARGE_POSITIVE_FLOAT, dt2,t=0;
51 err_t errorCode=SOLVER_NO_ERROR;
52 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
53 const dim_t blockSize=Paso_TransportProblem_getBlockSize(fctp);
54
55
56 if (dt<=0.) {
57 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");
58 }
59 if (blockSize>1) {
60 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: block size >0 is not supported.");
61 }
62 if (Paso_noError()) {
63 dt_max=Paso_TransportProblem_getSafeTimeStepSize(fctp);
64 /*
65 * allocate memory
66 *
67 */
68 fctfunction=Paso_FCTSolver_Function_alloc(fctp,options);
69 rsolver=Paso_ReactiveSolver_alloc(fctp);
70 u_save=TMPMEMALLOC(n,double);
71 Paso_checkPtr(u_save);
72 }
73 if (Paso_noError()) {
74 /*
75 * let the show begin!!!!
76 *
77 */
78 i_substeps=0;
79
80 /* while(i_substeps<n_substeps && Paso_noError()) */
81 if (dt_max < LARGE_POSITIVE_FLOAT) {
82 dt2=MIN(dt_max,dt);
83 } else {
84 dt2=dt;
85 }
86 Failed=0;
87
88 Paso_Copy(n,u,u0);
89 Paso_Copy(n,u_save, u);
90 while( (t<dt*(1.-sqrt(EPSILON))) && Paso_noError()) {
91
92 if (options->verbose) printf("substep step %d at t=%e (step size= %e)\n",i_substeps+1,t+dt2,dt2);
93
94 /* updates u */
95 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt2/2,q, options, &pp); /* Mu_t=Du+q u(0)=u */
96 if (errorCode == NO_ERROR) {
97 errorCode=Paso_FCTSolver_solve(fctfunction, u,dt2,options, &pp); /* Mv_t=Lv v(0)=u(dt/2) */
98 }
99 if (errorCode == NO_ERROR) {
100 errorCode=Paso_ReactiveSolver_solve(rsolver,fctp,u,dt2/2,q, options, &pp); /* Mu_t=Du+q u(dt/2)=v(dt/2) */
101 }
102 /*
103 * error handeling:
104 */
105 if (errorCode == SOLVER_NO_ERROR) {
106 Failed=0;
107 i_substeps++;
108 t+=dt2;
109 if (fctp->dt_max < LARGE_POSITIVE_FLOAT) {
110 dt2=MIN3(dt_max,dt2*1.1,dt-t);
111 } else {
112 dt2=MIN(dt2*1.1,dt-t);
113 }
114 Paso_Copy(n,u_save, u);
115 } else if ( (errorCode == SOLVER_MAXITER_REACHED) || (errorCode == SOLVER_DIVERGENCE) ) {
116 /* if FAILURES_MAX failures in a row: give up */
117 if (Failed > FAILURES_MAX) {
118 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: no convergence after time step reduction.");
119 } else {
120 if (options->verbose) printf("Paso_TransportProblem_solve: no convergence. Time step size is reduced.\n");
121 dt2*=0.5;
122 Failed++;
123 Paso_Copy(n,u, u_save);
124 }
125 } else if (errorCode == SOLVER_INPUT_ERROR) {
126 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: input error for solver.");
127 } else if (errorCode == SOLVER_MEMORY_ERROR) {
128 Paso_setError(MEMORY_ERROR,"Paso_TransportProblem_solve: memory allication failed.");
129 } else if (errorCode == SOLVER_BREAKDOWN) {
130 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: solver break down.");
131 } else if (errorCode == SOLVER_NEGATIVE_NORM_ERROR) {
132 Paso_setError(VALUE_ERROR,"Paso_TransportProblem_solve: negative norm.");
133 } else {
134 Paso_setError(SYSTEM_ERROR,"Paso_TransportProblem_solve: general error.");
135 }
136 }
137 }
138 /*
139 * clean-up:
140 *
141 */
142 Paso_FCTSolver_Function_free(fctfunction);
143 Paso_ReactiveSolver_free(rsolver);
144 TMPMEMFREE(u_save);
145
146 }
147 double Paso_TransportProblem_getSafeTimeStepSize(Paso_TransportProblem* fctp)
148 {
149 double dt_max, dt1, dt2;
150 if ( ! fctp->valid_matrices) {
151 /* set row-sum of mass_matrix */
152 Paso_SystemMatrix_rowSum(fctp->mass_matrix,fctp->lumped_mass_matrix);
153 /* split off row-sum from transport_matrix */
154 Paso_SystemMatrix_makeZeroRowSums(fctp->transport_matrix,fctp->reactive_matrix);
155 /* get a copy of the main diagonal of the mass matrix */
156 Paso_SystemMatrix_copyFromMainDiagonal(fctp->mass_matrix,fctp->main_diagonal_mass_matrix);
157
158
159 /* XXXXXXXXX hier geht es weiter XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX */
160 dt1=Paso_ReactiveSolver_getSafeTimeStepSize(fctp);
161 dt2=Paso_FCTSolver_getSafeTimeStepSize(fctp);
162 dt_max=MIN(dt1,dt2);
163
164 if (dt_max <= 0.) {
165 Paso_setError(TYPE_ERROR,"Paso_TransportProblem_solve: dt must be positive.");
166 }
167 fctp->dt_max=dt_max;
168 fctp->valid_matrices=Paso_noError();
169 }
170 return fctp->dt_max;
171 }

  ViewVC Help
Powered by ViewVC 1.1.26