source: trunk/sample/miscellaneous/Discrete/aborber_discrete.mso @ 978

Last change on this file since 978 was 978, checked in by Argimiro Resende Secchi, 6 years ago

examples of discrete polynomial approximation

File size: 11.1 KB
Line 
1#*---------------------------------------------------------------------
2* EMSO Model Library (EML) Copyright (C) 2004 - 2016 ALSOC.
3*
4* This LIBRARY is free software; you can distribute it and/or modify
5* it under the therms of the ALSOC FREE LICENSE as available at
6* http://www.enq.ufrgs.br/alsoc.
7*
8* EMSO Copyright (C) 2004 - 2016 ALSOC, original code
9* from http://www.rps.eng.br Copyright (C) 2002-2004.
10* All rights reserved.
11*
12* EMSO is distributed under the therms of the ALSOC LICENSE as
13* available at http://www.enq.ufrgs.br/alsoc.
14*
15*----------------------------------------------------------------------
16* Author: Argimiro R. Secchi
17* COPPE/UFRJ - Group of Modeling, Simulation, Control,
18*                               and Optimization of Processes
19* $Id$
20*--------------------------------------------------------------------*#
21
22FlowSheet Coluna_Abs
23       
24        PARAMETERS
25        Disc as Plugin(Type="Discrete", Boundary="BOTH",
26                                        InternalPoints=12, Lower=1, Upper=32,
27                                        Normal=32);
28
29        n  as Integer;
30        np as Integer;
31        Ns as Integer;
32        N1 as Integer;
33        n2 as Integer;
34        Am(Ns) as Real(Brief="recebe valores de do vetor Aplus");
35        An(Ns) as Real(Brief="recebe valores de do vetor Aminus");
36        Bm(N1) as Real(Brief="recebe valores de do vetor Aplus");
37        Bn(N1) as Real(Brief="recebe valores de do vetor Aminus");
38        r(np) as Real(Brief="raizes de hahn");
39        w(np) as Real(Brief="quadrature weights");
40        Norm  as Real(Brief="scaling");
41        V(n2) as Real(Brief="V matrix");
42       
43        Nest as Integer(Brief="Número de estágios");
44        alfa as Real;
45        beta as Real;
46        m        as Real;
47        h        as Real;
48        f        as Real;
49        B(N1) as Real;
50        A(Ns) as Real;
51       
52        SET
53        np= Disc.NodalPoints;
54        n = Disc.InternalPoints;
55        n2 = 2*n;
56        Ns= np * np;
57        N1= n*np;
58        Am = Disc.Aplus;       
59        An = Disc.Aminus;
60        r = Disc.Roots;
61        w = Disc.Weights;
62        Norm = Disc.Scaling;
63        V = Disc.Vmatrix;
64        Bm = Disc.Bplus;       
65        Bn = Disc.Bminus;
66
67        Nest = Disc.UpperBound-Disc.LowerBound+1;
68        alfa = 0.75;
69        m = 2;
70        h = 6;
71        f = 0;
72        beta = alfa + m;
73        B = alfa * Bn + m * Bm;
74        A = alfa * An + m * Am;
75
76        VARIABLES
77       
78        rs(np) as Real(Brief="desnormal.");
79        X(Nest) as Real(Brief="Perfil completo");
80        Xred(np) as Real(Brief="Perfil momentos modificado");
81        Xmom(np) as Real(Brief="Perfil momentos nos pontos de colocação");
82        Xcol(np) as Real(Brief="Perfil colocação nos pontos de colocação");
83
84        Xmint(Nest) as Real(Brief="Perfil momentos interpolado");
85        Xcint(Nest) as Real(Brief="Perfil colocação interpolado");
86        Y(Nest) as Real;
87        Ymom(Nest) as Real;
88        Ycol(Nest) as Real;
89
90        Em as Real(Brief="Integral do erro quadrático médio dos momentos");
91        Ec as Real(Brief="Integral do erro quadrático médio da colocação");
92
93        Emi(Nest) as Real(Brief="Integral do erro dos momentos");
94        Eci(Nest) as Real(Brief="Integral do erro da colocação");
95
96        EQUATIONS
97       
98        rs = r * Norm;
99
100# Completa
101
102        X(1) = f;
103        X(Nest) = h/m;
104
105        for i in [2:Nest-1] do
106           diff(X(i))*'s' = alfa*X(i-1)+m*X(i+1)-beta*X(i);
107        end
108
109        Y = m * X;
110
111# Momentos
112
113        Xmom(1) = f;
114        Xmom(np) = h/m;
115        Xred(1) = Xmom(1);
116        Xred(np) = Xmom(np);
117        Xred(2:n+1) = Xmom(2:n+1)+V(1:n)*Xmom(1)+V(n+1:2*n)*Xmom(np);
118       
119        for i in [2:n+1] do
120           diff(Xred(i))*'s' = sum(B([1:np]+(i-2)*np)*Xmom)-beta*Xred(i);
121        end
122
123# Colocação ortogonal
124
125        Xcol(1) = f;
126        Xcol(np) = h/m;
127
128        for i in [2:n+1] do
129           diff(Xcol(i))*'s' = sum(A([1:np]+(i-1)*np)*Xcol)-beta*Xcol(i);
130        end
131
132# Interpolations
133        for i in [1:Nest] do
134                Xmint(i) = Disc.Interpol(Xmom,i);
135                Xcint(i) = Disc.Interpol(Xcol,i);
136        end
137
138        Ymom = m * Xmint;
139        Ycol = m * Xcint;
140
141# Cálculo dos erros
142
143        diff(Em)*'s' = sqrt(sum((X - Xmint)^2))/Nest;
144        diff(Ec)*'s' = sqrt(sum((X - Xcint)^2))/Nest;
145
146        diff(Emi)*'s' = abs(X - Xmint);
147        diff(Eci)*'s' = abs(X - Xcint);
148
149        INITIAL
150        X(2:Nest-1) = 0;
151        Xred(2:n+1) = 0+V(1:n)*Xmom(1)+V(n+1:2*n)*Ymom(1);
152        Xcol(2:n+1) = 0;
153        Em = 0;
154        Ec = 0;
155        Emi = 0;
156        Eci = 0;
157
158        OPTIONS
159        Dynamic         = true; 
160        TimeStep        = 0.2;
161        TimeEnd         = 30; # colocar igual ao número de estágios internos
162        DAESolver(File="dasslc", RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
163        NLASolver(File="sundials", RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
164end
165
166FlowSheet Coluna_Abs_nl
167       
168        PARAMETERS
169        Disc as Plugin(Type="Discrete", Boundary="BOTH",
170                                        InternalPoints=5, Lower=1, Upper=32,
171                                        Normal=32);
172
173        n  as Integer;
174        np as Integer;
175        Ns as Integer;
176        N1 as Integer;
177        n2 as Integer;
178        Am(Ns) as Real(Brief="recebe valores de do vetor Aplus");
179        An(Ns) as Real(Brief="recebe valores de do vetor Aminus");
180        Bm(N1) as Real(Brief="recebe valores de do vetor Aplus");
181        Bn(N1) as Real(Brief="recebe valores de do vetor Aminus");
182        r(np) as Real(Brief="raizes de hahn");
183        w(np) as Real(Brief="quadrature weights");
184        Norm  as Real(Brief="scaling");
185        V(n2) as Real(Brief="V matrix");
186       
187        Nest as Integer(Brief="Número de estágios");
188        alfa as Real;
189        beta as Real;
190        nu   as Real;
191        h        as Real;
192        f        as Real;
193       
194        SET
195        np= Disc.NodalPoints;
196        n = Disc.InternalPoints;
197        n2 = 2*n;
198        Ns= np * np;
199        N1= n*np;
200        Am = Disc.Aplus;       
201        An = Disc.Aminus;
202        r = Disc.Roots;
203        w = Disc.Weights;
204        Norm = Disc.Scaling;
205        V = Disc.Vmatrix;
206        Bm = Disc.Bplus;       
207        Bn = Disc.Bminus;
208
209        Nest = Disc.UpperBound-Disc.LowerBound+1;
210        beta = 0.75;
211        alfa = 1.2;
212        nu = alfa - 1;
213        h = 0.9; # 1.2
214        f = 0;
215
216        VARIABLES
217       
218        rs(np) as Real(Brief="desnormal.");
219        X(Nest) as Real(Brief="Perfil completo");
220        Xred(np) as Real(Brief="Perfil momentos modificado");
221        Yred(np) as Real(Brief="Perfil momentos modificado");
222        Xmom(np) as Real(Brief="Perfil momentos nos pontos de colocação");
223        Ymom(np) as Real(Brief="Perfil momentos nos pontos de colocação");
224        Xcol(np) as Real(Brief="Perfil colocação nos pontos de colocação");
225        Ycol(np) as Real(Brief="Perfil colocação nos pontos de colocação");
226
227        Xmint(Nest) as Real(Brief="Perfil momentos interpolado");
228        Xcint(Nest) as Real(Brief="Perfil colocação interpolado");
229        Y(Nest) as Real;
230        Ymint(Nest) as Real;
231        Ycint(Nest) as Real;
232
233        Em as Real(Brief="Integral do erro quadrático médio dos momentos");
234        Ec as Real(Brief="Integral do erro quadrático médio da colocação");
235
236        Emi(Nest) as Real(Brief="Integral do erro dos momentos");
237        Eci(Nest) as Real(Brief="Integral do erro da colocação");
238
239        EQUATIONS
240       
241        rs = r * Norm;
242
243# Completa
244
245        X(1) = f;
246        Y(Nest) = h;
247
248        for i in [2:Nest-1] do
249           diff(X(i))*'s' = beta*X(i-1)+Y(i+1)-beta*X(i)-Y(i);
250        end
251
252        Y = alfa * X / (1 + nu * X);
253
254# Momentos
255
256        Xmom(1) = f;
257        Ymom(np) = h;
258        Xred(1) = Xmom(1);
259        Xred(np) = Xmom(np);
260        Xred(2:n+1) = Xmom(2:n+1)+V(1:n)*Xmom(1)+V(n+1:2*n)*Xmom(np);
261        Yred(1) = Ymom(1);
262        Yred(np) = Ymom(np);
263        Yred(2:n+1) = Ymom(2:n+1)+V(1:n)*Ymom(1)+V(n+1:2*n)*Ymom(np);
264       
265        for i in [2:n+1] do
266           diff(Xred(i))*'s' = sum(Bm([1:np]+(i-2)*np)*Ymom)+
267                                                beta*sum(Bn([1:np]+(i-2)*np)*Xmom)-
268                                                beta*Xred(i)-Yred(i);
269        end
270
271        Ymom = alfa * Xmom / (1 + nu * Xmom);
272
273# Colocação ortogonal
274
275        Xcol(1) = f;
276        Ycol(np) = h;
277
278        for i in [2:n+1] do
279           diff(Xcol(i))*'s' = sum(Am([1:np]+(i-1)*np)*Ycol)+
280                                                beta*sum(An([1:np]+(i-1)*np)*Xcol)-
281                                                beta*Xcol(i)-Ycol(i);
282        end
283
284        Ycol = alfa * Xcol / (1 + nu * Xcol);
285
286# Interpolations
287        for i in [1:Nest] do
288                Xmint(i) = Disc.Interpol(Xmom,i);
289                Xcint(i) = Disc.Interpol(Xcol,i);
290                Ymint(i) = Disc.Interpol(Ymom,i);
291                Ycint(i) = Disc.Interpol(Ycol,i);
292        end
293
294# Cálculo dos erros
295
296        diff(Em)*'s' = sqrt(sum((X - Xmint)^2))/Nest;
297        diff(Ec)*'s' = sqrt(sum((X - Xcint)^2))/Nest;
298
299        diff(Emi)*'s' = abs(X - Xmint);
300        diff(Eci)*'s' = abs(X - Xcint);
301
302        INITIAL
303        X(2:Nest-1) = 0;
304        Xred(2:n+1) = 0+V(1:n)*Xmom(1)+V(n+1:2*n)*Xmom(np);
305        Xcol(2:n+1) = 0;
306        Em = 0;
307        Ec = 0;
308        Emi = 0;
309        Eci = 0;
310
311        OPTIONS
312        Dynamic         = true; 
313        TimeStep        = 0.2;
314        TimeEnd         = 30; # colocar igual ao número de estágios internos
315#       DAESolver(File="dasslc", RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
316#       NLASolver(File="sundials", RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);   
317end
318
319FlowSheet Coluna_Abs_nl_ss
320       
321        PARAMETERS
322        Disc as Plugin(Type="Discrete", Boundary="BOTH",
323                                        InternalPoints=5, Lower=1, Upper=32,
324                                        Normal=32);
325
326        n  as Integer;
327        np as Integer;
328        Ns as Integer;
329        N1 as Integer;
330        n2 as Integer;
331        Am(Ns) as Real(Brief="recebe valores de do vetor Aplus");
332        An(Ns) as Real(Brief="recebe valores de do vetor Aminus");
333        Bm(N1) as Real(Brief="recebe valores de do vetor Aplus");
334        Bn(N1) as Real(Brief="recebe valores de do vetor Aminus");
335        r(np) as Real(Brief="raizes de hahn");
336        w(np) as Real(Brief="quadrature weights");
337        Norm  as Real(Brief="scaling");
338        V(n2) as Real(Brief="V matrix");
339       
340        Nest as Integer(Brief="Número de estágios");
341        alfa as Real;
342        beta as Real;
343        nu   as Real;
344        h        as Real;
345        f        as Real;
346       
347        SET
348        np= Disc.NodalPoints;
349        n = Disc.InternalPoints;
350        n2 = 2*n;
351        Ns= np * np;
352        N1= n*np;
353        Am = Disc.Aplus;       
354        An = Disc.Aminus;
355        r = Disc.Roots;
356        w = Disc.Weights;
357        Norm = Disc.Scaling;
358        V = Disc.Vmatrix;
359        Bm = Disc.Bplus;       
360        Bn = Disc.Bminus;
361
362        Nest = Disc.UpperBound-Disc.LowerBound+1;
363        beta = 1; #1.25;
364        alfa = 0.75; #1.2;
365        nu = alfa - 1;
366        h = 1; # 1.2
367        f = 0;
368
369        VARIABLES
370       
371        rs(np) as Real(Brief="desnormal.");
372        X(Nest) as Real(Brief="Perfil completo");
373        Xred(np) as Real(Brief="Perfil momentos modificado");
374        Yred(np) as Real(Brief="Perfil momentos modificado");
375        Xmom(np) as Real(Brief="Perfil momentos nos pontos de colocação");
376        Ymom(np) as Real(Brief="Perfil momentos nos pontos de colocação");
377        Xcol(np) as Real(Brief="Perfil colocação nos pontos de colocação");
378        Ycol(np) as Real(Brief="Perfil colocação nos pontos de colocação");
379
380        Xmint(Nest) as Real(Brief="Perfil momentos interpolado");
381        Xcint(Nest) as Real(Brief="Perfil colocação interpolado");
382        Y(Nest) as Real;
383        Ymint(Nest) as Real;
384        Ycint(Nest) as Real;
385
386        EQUATIONS
387       
388        rs = r * Norm;
389
390# Completa
391
392        X(1) = f;
393        Y(Nest) = h;
394
395        for i in [2:Nest-1] do
396           diff(X(i))*'s' = beta*X(i-1)+Y(i+1)-beta*X(i)-Y(i);
397        end
398
399        Y = alfa * X / (1 + nu * X);
400
401# Momentos
402
403        Xmom(1) = f;
404        Ymom(np) = h;
405        Xred(1) = Xmom(1);
406        Xred(np) = Xmom(np);
407        Xred(2:n+1) = Xmom(2:n+1)+V(1:n)*Xmom(1)+V(n+1:2*n)*Xmom(np);
408        Yred(1) = Ymom(1);
409        Yred(np) = Ymom(np);
410        Yred(2:n+1) = Ymom(2:n+1)+V(1:n)*Ymom(1)+V(n+1:2*n)*Ymom(np);
411       
412        for i in [2:n+1] do
413           diff(Xred(i))*'s' = sum(Bm([1:np]+(i-2)*np)*Ymom)+
414                                                beta*sum(Bn([1:np]+(i-2)*np)*Xmom)-
415                                                beta*Xred(i)-Yred(i);
416        end
417
418        Ymom = alfa * Xmom / (1 + nu * Xmom);
419
420# Colocação ortogonal
421
422        Xcol(1) = f;
423        Ycol(np) = h;
424
425        for i in [2:n+1] do
426           diff(Xcol(i))*'s' = sum(Am([1:np]+(i-1)*np)*Ycol)+
427                                                beta*sum(An([1:np]+(i-1)*np)*Xcol)-
428                                                beta*Xcol(i)-Ycol(i);
429        end
430
431        Ycol = alfa * Xcol / (1 + nu * Xcol);
432
433# Interpolations
434        for i in [1:Nest] do
435                Xmint(i) = Disc.Interpol(Xmom,i);
436                Xcint(i) = Disc.Interpol(Xcol,i);
437                Ymint(i) = Disc.Interpol(Ymom,i);
438                Ycint(i) = Disc.Interpol(Ycol,i);
439        end
440
441        INITIAL
442        X(2:Nest-1) = 0;
443        Xred(2:n+1) = 0+V(1:n)*Xmom(1)+V(n+1:2*n)*Xmom(np);
444        Xcol(2:n+1) = 0;
445
446        OPTIONS
447        Dynamic         = false; 
448        TimeStep        = 0.2;
449        TimeEnd         = 30; # colocar igual ao número de estágios internos
450#       DAESolver(File="dasslc", RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);
451#       NLASolver(File="sundials", RelativeAccuracy = 1e-6,AbsoluteAccuracy = 1e-8);   
452end
Note: See TracBrowser for help on using the repository browser.