JP2005267418A - Numerical calculation method, numerical calculation device and numerical calculation program - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To enhance convergent property in a numerical calculation using a residual removal method adapting AMG (algebraic multi grid) method for an internal solver. <P>SOLUTION: A variable separation type AMG method is used. As an AMG cycle, v-cycle is used. In the v-cycle, in case of descending the grid level, the repetition frequency n of moderation is set relatively small when the level is shallow, and set relatively large when the level is deep. On the other hand, in case of ascending the grid level, the repetition frequency n of moderation is set to a predetermined value n<SB>u</SB>. According to this, a result close to the V-cycle can be obtained with a minimized calculation quantity. <P>COPYRIGHT: (C)2005,JPO&NCIPI

Description

  The present invention relates to a numerical calculation technique using a computer.

  Conventionally, successive approximation algorithms have been used in control and numerical calculations for physical quantities such as pressure and temperature. The successive approximation algorithm is an algorithm for obtaining an optimal solution by numerical calculation while repeatedly performing successive approximation based on various optimal control theories.

  A residual ablation method is known as a highly convergent successive approximation algorithm that can obtain a solution for each target parameter at high speed and with high accuracy. And in patent document 1, the technique which implement | achieves the further improvement of convergence in the residual cutting method is proposed.

  As one of acceleration methods of relaxation methods (Gauss Seidel method, SOR method, etc.), there is a multi grid (MG) method. In general, the relaxation method can quickly remove the high-frequency component of the residual, but requires a very large number of iterations to attenuate the low-frequency component. In the MG method, a plurality of lattices having different lattice intervals (roughness) are prepared, and by using a relaxation method on each lattice, a high-frequency component of a residual in each lattice is removed, and a solution is obtained at high speed. The low frequency component of the residual in a certain grid becomes a high frequency component in a coarser grid.

As a derivation method of the MG method, there is an algebraic multigrid (AMG) method. In order to use the above-described MG method (sometimes called a geometric multigrid method in order to distinguish it from the AMG method), it is necessary to prepare a lattice system manually in advance. On the other hand, when the AMG method is used, it is not necessary to provide a lattice system, and it is only necessary to provide simultaneous linear equations. In the AMG method, a virtual lattice is created from a given coefficient matrix of simultaneous linear equations, and a sparse grid is also created only from the information of the coefficient matrix.
JP 2001-134304 A Andrew J. Cleary, et al. "ROBUSTNESS AND SCALABILITY OF ALGEBRAIC MULTIGRID", SIAM J. SCI. COMPUT., 2000, Vol.21, No.5, pp.1886-1908

  The inventors of the present application are considering using the above-described AMG method as an internal solver of the residual cutting method.

  On the other hand, in simultaneous linear equations obtained from structural analysis and fluid analysis, unknown variables may include a plurality of physical quantities. In such a case, since the relationship between physical quantities is added to the relationship between lattices, the coefficient matrix of the simultaneous linear equations to be solved contains more non-zero elements than in the case of a single physical quantity. For this reason, problems are often difficult to solve (stiff).

  For such a stiff problem, simply applying the residual cutting method adopting the AMG method to the internal solver as described above, a convergent solution cannot be obtained, or a long calculation time is required until a converged solution is obtained. May be required.

  In view of the above problems, the present invention has an object to improve convergence in numerical calculation using a successive approximation algorithm by a computer, in particular, numerical calculation using a residual cutting method employing the AMG method as an internal solver. And

In order to solve the above-mentioned problem, the solution means taken by the first invention is that a coefficient matrix (N rows and N columns, N is a positive integer) obtained by discretizing a partial differential equation to be satisfied by the physical quantity U is A Assuming that the inhomogeneous term (source term) is f, a method, apparatus, and program for solving a numerical value of the physical quantity U by solving A · U = f using a computer, an initial value U of the physical quantity U ⁰ is set, 0 is given as the initial value to the number of iterations m, 0 is given as the initial value of the perturbation amount φ, (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r, and the number of iterations m while incrementing the predicted approximate value [psi ^m of a · φ = r ^m, and the minimum a first step of obtaining by iteration by the first calculation unit having an internal solver, the L2 norm of the residual r ^{m + 1} The corrected approximate value φ ^m to be Calculated from the predicted approximate value [psi ^m Te, given as an approximate solution ^{U m + 1 (U m +} φ m), a second step of providing a residual r ^{m + 1} a ^{(r m -A · φ m)} , The approximate solution U ^m is repeatedly executed until it converges, and the first step uses the AMG method for the internal solver, and as the AMG cycle, the number of relaxation iterations n determines the grid level. When descending, the level is relatively small when the level is relatively shallow, and relatively large when the level is relatively deep. On the other hand, when increasing the grid level, ν− is set to a predetermined value. A cycle is used.

In the ν-cycle according to the first aspect of the present invention, n = max (k−n _b , 0) is applied to the relaxation performed before the number of iterations n of the relaxation method is “limited” from 1) level k to level k + 1. 2) For relaxation performed before “interpolating” from level k + 1 to level k, n = _nu is preferably set (n _b and n _u are natural numbers).

The solution provided by the second invention is that the coefficient matrix (N rows and N columns, N is a positive integer) when discretizing the partial differential equation to be satisfied by the physical quantity U is A, an inhomogeneous term (source term) ) Is a method, apparatus, and program for solving a numerical value of the physical quantity U by solving A · U = f using a computer, and setting an initial value U ⁰ of the physical quantity U, 0 is given as an initial value to the number of iterations m, 0 is given as the initial value of the perturbation amount φ, (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r, and the iteration number m is incremented, the predicted approximate value [psi ^m of a · φ = r ^m, a first step of obtaining by iteration by the first calculation unit having an internal solver, modified approximation which minimizes the L2 norm of the residual r ^{m + 1} phi the ^m, from the predicted approximate value [psi ^m by the optimization routine by the second arithmetic unit Because, given as an approximate solution ^{U m + 1 (U m +} φ m), a second step of providing a residual r ^{m + 1} a ^{(r m -A · φ m)} , the approximate solution U ^m is converged When the residual cutting rate becomes smaller than a predetermined value, the modified approximate solution φ ^m is discarded, and a reset process that involves changing the algorithm of the internal solver used in the first step is performed. Is what you do.

  In the second invention, the first step uses an AMG method as the internal solver. In the reset process, as the algorithm of the internal solver is changed, the number of iterations n of relaxation in the AMG cycle is n. It is preferable to change the setting.

  According to the first invention, the ν-cycle is used as the AMG cycle in the internal solver of the residual cutting method. In the ν-cycle, when descending the grid level, the number of iterations of relaxation is set to a low level at a shallow level where the influence of the relaxation performed before “restriction” is considered to be small, and a large number is set to a deep level. On the other hand, when increasing the grid level, the number of iterations of relaxation after “interpolation” is set to a predetermined value as in the V-cycle. Therefore, it is possible to obtain almost the same performance as the V-cycle with a smaller amount of calculation.

  According to the second invention, when the residual cutting rate becomes smaller than a predetermined value, the reset process is performed and the algorithm of the internal solver is changed. As a result, it is possible to increase the possibility of finding a solution correction direction that is different from the previous one, since the ability to find the solution correction direction is lost. Therefore, the convergence can be further improved.

  Hereinafter, embodiments of the present invention will be described in detail with reference to the drawings.

  The numerical calculation method according to the present embodiment is based on the residual cutting method using the AMG method for the internal solver.

<AMG method>
First, an outline of the AMG method will be described.

  The processing of the AMG method is roughly divided into two types: AMG setup and AMG cycle. In the AMG setup, generation of sparse grids, generation of “restrictions” and “interpolation” functions between sparse grids, and generation of coefficient matrices in sparse grids are performed. Since these details are described in detail in Non-Patent Document 1 described above, the description thereof is omitted here.

(AMG cycle)
The AMG cycle is a phase in which simultaneous linear equations are solved using a sparse grid generated in an AMG setup, a “restriction” “interpolation” function, and a coefficient matrix. The simplest cycle is a case where there are only two grids, and the following procedures A to E result in one cycle.

A. Apply n-time relaxation method (Gauss-Seidel method, etc.) to A ⁰ U ⁰ = f ⁰ (relaxation before restriction)
B. Using the limiting function I ₀ ¹ , f ¹ = I ₀ ¹ (f ⁰ −A ⁰ U ⁰ ) is calculated.

C. Solve A ¹ U ¹ = f ^1, obtaining the U ^1.

D. Using an interpolation function I ₁ ^0, corrects _{^{U 0 ← U 0 + I 1}} 0 U 1 of the approximate solution.

E. Apply the relaxation method n times to A ⁰ U ⁰ = f ⁰ (relaxation after interpolation).

  The superscript here indicates the grid level. An expression with a subscript of 0, that is, an expression at grid level 0 corresponds to the above-described expression A · U = f without a subscript. Note that, although confusing, the superscripts other than the AMG cycle indicate the number of iterations.

  The cycle of the above procedures A to E is expressed as shown in FIG. In FIG. 1, the vertical direction is the grid level, and the horizontal direction is the processing order. Procedures A and E are processing at the same grid level 0, and procedure C is processing at grid level 1. The procedure B is a process of transitioning the grid level from the level 0 to the level 1 by the “restriction” function, and the procedure D is a process of transitioning the grid level from the level 1 to the level 0 by the “interpolation” function. As shown in FIG. 1, such processing is called “V-cycle” because the shape represented by dots and lines is similar to the letter “V”.

  The V-cycle when the grid has three or more stages is the same as the case of two stages. That is, the grid level is lowered by repeating the procedures A and B, the solution is obtained at the deepest grid (procedure C), the grid level is raised by repeating the procedures D and E, and one cycle is completed.

  In addition to V-cycles, W-cycles and sawtooth-cycles have been proposed. FIG. 2 is a conceptual diagram showing a W-cycle (a) and a sawtooth-cycle (b) when the grid has three stages.

(Variable separation type AMG method)
As described above, in the AMG method, sparse grid generation (lattice thinning) is performed based on information called a coefficient matrix. When there is only one kind of unknown in the differential equation that is the basis of the simultaneous linear equations to be solved, that is, when there is only one unknown on one grid, it is often possible to infer grid information from a coefficient matrix. . However, when dealing with multivariate problems, it is generally difficult to read grid information from only the coefficient matrix. A technique devised to avoid this difficulty is the variable separation AMG method.

  In the variable separation type AMG method, in order to generate a sparse grid, in addition to the coefficient matrix, information on the type of unknown amount is used. The specific procedure is as follows.

  A. In the dense grid, the unknown quantities are separated for each type, and the simultaneous linear equations are constructed for each type of unknown quantities by ignoring the terms representing the interactions between the different kinds of unknown quantities in the coefficient matrix.

  B. Using a general AMG setup algorithm, a sparse grid and a “restriction” and “interpolation” function between sparse and dense grids are generated for each type of unknown quantity.

  C. Constructing a coefficient matrix of simultaneous linear equations on a sparse grid using the sparse grid generated in step B, the "constraint" "interpolation" function, and the original coefficient matrix that does not ignore the terms representing the interaction between unknowns To do.

  FIG. 3 is a flowchart showing a numerical calculation method according to this embodiment. FIG. 4 is a block diagram showing the configuration of the numerical calculation apparatus 1 according to this embodiment. In the present embodiment, the above-described variable separation type AMG method is used as an internal solver of the residual cutting method for obtaining the solution correction vector Ψ. This numerical calculation method can be implemented by causing a computer to execute a program for realizing the method, which is recorded on a recording medium such as the CD-ROM 2. In addition, it can be easily realized by an apparatus including a computer that executes a program for realizing the method.

  The numerical calculation method according to the present embodiment includes some features that are not present in the past. One of them is the “ν-AMG cycle”, and the other is “change of the AMG cycle when a reset occurs”.

<Ν-AMG cycle>
In this embodiment, it is assumed that only one cycle as described below is used in one process in the internal solver.

Assume that the level of the grid to which the unknown variable belongs initially is level 0, and the grid number increases as level 1, level 2,. When n _b and n _u are natural numbers and the number of iterations of the relaxation method is n,
1) For mitigation before “restricting” from level k to level k + 1:
n = max (k−n _b , 0)
2) For mitigation before “interpolating” from level k + 1 to level k:
n = n _u
And

FIG. 5 shows an example of the AMG cycle used in this embodiment. In FIG. 5, a 6-stage grid is used, and n _b = 1. In the restriction from level 0 to level 1 and the restriction from level 1 to level 2, the number of relaxation iterations n = 0. The present inventor has named such an AMG cycle as a “ν-cycle”.

  The advantages of the ν-cycle devised by the present inventor will be described.

  The role of the internal solver in the residual cutting method is to find a direction independent of the solution correction vector sequence of the past iterations, and the residual equation does not need to be accurately solved by the internal solver. In extreme terms, the residual solver internal solver may not have the ability to solve simultaneous linear equations alone.

  When the AMG method was used for the internal solver, the following was found from the results of numerical experiments. When a cycle with more procedures than the V-cycle (such as a W-cycle) is used, there is no significant change in the convergence rate per iteration of the residual cutting method regardless of which cycle is used. On the other hand, when the sawtooth cycle is used, the convergence rate is significantly reduced. In the V-cycle, it is possible to find a direction independent of the solution correction vector sequence in the past, but it is often difficult to find in the sawtooth cycle. The ν-cycle is closer to the sawtooth-cycle than the V-cycle, but requires less computation than the sawtooth-cycle. This is because in the sawtooth cycle, the relaxation that occurs before “limiting” from the closest grid to the coarser grid is not performed in the v-cycle.

  And, for the following reason, it is considered that a good result can be obtained when the ν-cycle having a smaller calculation amount than the sawtooth-cycle is used. That is, since the problem calculated by the internal solver of the residual cutting method often uses a zero vector as the initial value, it is considered that the effect of relaxation performed before “restriction” is small when the grid level is shallow. In the v-cycle, at shallow grid levels where the effect of relaxation is considered to be small, relaxation is not performed before “restriction” or the number of relaxation iterations is reduced. On the other hand, at the deep grid level, the number of relaxation iterations performed before the “limit” is increased. In addition, the relaxation after “interpolation” when raising the grid level is performed by setting the number of iterations to a fixed number, as in the V-cycle. Therefore, the ν-cycle is considered to produce results close to the V-cycle.

<AMG cycle change at reset>
For the relative residual ε ^m (= ‖r ^m ‖ / ‖b‖) in the m-th iteration, the residual cutting rate δ ^m = (ε ^m −ε ^{m + 1} ) / ε ^m is larger than a predetermined value Δ. When it becomes smaller, the following processing is performed.

Process A The sequence φ ^l (l = m−L _max +1,..., M) of solution correction vectors in the past iterations up to the mth iteration is discarded.

  Process B Change the internal solver algorithm.

  Process A is a technique commonly known as iterative reset. It is known as a coping method when the iteration method falls into a local solution, and a method of resetting at a predetermined number of iterations and a method of resetting at every predetermined number of iterations have been reported.

  In the residual cutting method, an internal solver detects the direction (vector) in which the solution should be corrected, and based on the assumption that the approximate solution is represented by a linear combination of the solution correction vector and the solution correction vector of the past iteration. The coefficients of the linear combination are determined so that the relative residual in the next iteration is minimized. Here, the fact that the above-described residual cutting rate is extremely small is considered that the internal solver used does not have the ability to find a direction in which the solution should be corrected any more. For this reason, changing the algorithm of the internal solver increases the possibility of finding a different solution correction direction.

  In this embodiment, as an example of changing the algorithm, the setting of the number of relaxation iterations n in the AMG cycle is changed.

Specifically, for example, changes are made as follows. That is, the AMG algorithm using the ν-cycle as shown in FIG. 5 is used from the start of calculation until the first reset occurs. After the first reset, the number of iterations of relaxation performed before “limit” and after “interpolation” is all n = _nu . After the second and subsequent resets, a natural number _ns is added to the number of iterations n so far, and this is set as a new iteration number n. When n reaches a predetermined value N, all the processing is terminated as non-convergence.

  The advantage of the method of changing the AMG cycle at the time of reset devised by the present inventor is considered as follows. As discussed in the above-mentioned Patent Document 1, it is difficult and important to determine the number of relaxations in the internal solver in the residual cutting method. If the number of relaxations is reduced, the amount of calculation per iteration is reduced, but the convergence rate is deteriorated, so that there is a high possibility that the convergence solution cannot be reached. Conversely, when the number of relaxations is increased, the amount of calculation per iteration increases, but the convergence rate increases and a converged solution can be obtained stably.

  In the change of the AMG cycle described above, the calculation is initially performed using the ν-cycle with a small amount of calculation, and when the convergence rate (residual cutting rate) becomes small, the amount of calculation is large but the ability to find the solution correction direction is high. When the V-cycle is considered to be changed and the convergence rate is further reduced, the number of relaxations at each stage of the V-cycle is increased. As a result, a convergent solution can be obtained stably with as little calculation as possible.

  The numerical calculation method according to the present embodiment will be described with reference to the flowchart shown in FIG.

If the physical quantity to be obtained is U, the coefficient matrix when the partial differential equation to be satisfied by the physical quantity U is discretized is A (N rows and N columns), and the inhomogeneous term (source term) is f, the equation to be solved is ,
A ・ U = f (1)
It is represented by In general, this equation can be solved by a successive approximation method such as the SOR method or the ADI method.

The solution U∞ of Equation (1) is expressed as follows by the approximate solution U and the perturbation amount (difference from the true solution) φ.
U∞ = U + φ (2)
Then, the true solution U∞ of Equation (1) is obtained by obtaining the perturbation amount φ and correcting the approximate solution U.

Here, the residual r for the approximate solution U is defined as follows.
r = f−A · U (3)
From formulas (1) to (3),
A · (U + φ) = f
Ａ A ・ φ = f−A ・ U
= R
Therefore, in order to obtain the perturbation amount φ, the following equation should be solved.
A · φ = r (4)

  In the present embodiment, instead of obtaining the convergence solution of Equation (4), the predicted approximate value Ψ is obtained by the iteration of the smallest unit with the steepest convergence gradient, and the obtained predicted approximate value Ψ is supplied to the optimization control routine. To do. The calculation of the predicted approximate value Ψ is repeatedly executed until the execution result of the optimization control routine satisfies a predetermined condition. Details of the processing are described in Patent Document 1.

First, in step S1, an initial value U ⁰ of the target physical quantity U is set. In step S2, AMG setup is performed. That is, generation of sparse grids, generation of “restrictions” and “interpolation” functions between sparse and dense grids, and generation of coefficient matrices in sparse grids are performed.

Next, in step S3, 0 is given as the initial value to the number of iterations m, 0 is given as the initial value of the perturbation amount φ, and (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r. Also, ‖r‖ / ‖f‖ is given as the initial value ε ⁰ of the relative residual ε. Further, regarding the number of relaxation iterations in the AMG cycle, n = 0 and n _u and n _s are set to predetermined values. In addition, a predetermined value Δ for determining reset is set.

Thereafter, steps S4 to S8 are repeatedly executed while incrementing the number of repetitions m until the approximate solution U ^m converges (S9, S11). When the residual cutting rate δ ^m becomes smaller than the predetermined value Δ, the modified approximate solution φ ^m is discarded and the algorithm of the internal solver is changed. Here, it is assumed that the setting of the number n of iterations of relaxation in the AMG cycle is changed.

In steps S4 to S7 (first step), the first arithmetic unit 10 having the internal solver 11 employing the AMG method,
A ・ φ = r ^m
A predicted approximate value Ψ ^{m of} is obtained. Since n = 0 is set from the beginning of the calculation until the first reset is performed (No in step S4), the above-described ν-cycle is executed as an AMG cycle in step S5. That is, the number of iterations n of the relaxation method is
1) For mitigation before “restricting” from level k to level k + 1:
n = max (k−n _b , 0)
2) For mitigation before “interpolating” from level k + 1 to level k:
n = n _u
( _Nb is a natural number).

When the predicted approximate value Ψ ^m is obtained in steps S4 to S7, in step S8, using the predicted approximate value Ψ ^m , L2 of the next residual r ^{m + 1} is obtained by the optimization routine in the second arithmetic unit 20. An optimized modified approximate value φ ^m is obtained so that the norm is minimized. The processing in the optimization routine such as calculation of the residual minimization coefficient α _l (l = 1,..., L) is described in detail in Patent Document 1, and is omitted here.

Then, the approximate value U ^{m + 1} and the residual r ^{m + 1} in the (m + 1) th iteration are obtained from the optimized modified approximate value φ ^m by the following equation.

U ^{m + 1} = U ^m + φ ^m
r ^{m + 1} = r ^m −A · φ ^m
Further, a relative residual ε ^{m + 1} and a residual cutting rate δ ^m are also obtained.

When the approximate solution U ^m has not yet converged (No in step S9), the residual cutting rate δ ^m is compared with a predetermined value Δ (step S10). If the residual excision rate δ ^m exceeds the predetermined value Δ (Yes in S10), the repeat count m is incremented (S11), and the process returns to step S4.

On the other hand, when the residual cutting rate δ ^m falls below the predetermined value Δ (No in S10), reset processing is performed to change the algorithm of the internal solver.

At the time of the first reset, since n = 0 (No in step S12), n = _nu is set. In addition, the iterative past solution correction vector φ ^l (l = m−L _max +1,..., M) is set to zero (step S13). Thereby, in subsequent iterations, the process proceeds to step S6 in the branch of step S4, and the V-cycle is executed as an AMG cycle in the internal solver. The number of relaxation iterations at each stage of the V-cycle is n = _nu .

In the second and subsequent resets, since n> 0 (Yes in step S12), n _s is added to n _u (step S14). Then, n = _nu . In addition, the iterative past solution correction vector φ ^l (l = m−L _max +1,..., M) is set to zero (step S13). Thereby, in subsequent iterations, the AMG method by the V-cycle in which the number of iterations of relaxation is increased by n _s in the internal solver is executed.

  Here, the internal solver algorithm at the time of resetting is changed by changing the setting of the number of iterations of relaxation in the AMG cycle. However, the present invention is not limited to this, and is completely different from a certain algorithm, for example. It is possible to change to this algorithm. Further, the method of changing the number of iterations of relaxation is not limited to the one shown here. For example, a V-cycle may be adopted from the beginning to increase the number of relaxations for each reset.

  Further, even when the reset process is not performed, by adopting the ν-cycle, it is possible to obtain almost the same performance as the V-cycle with a small amount of calculation.

<Numerical experiment>
Using the numerical calculation method shown in this embodiment, a simultaneous linear equation obtained by fluid analysis was actually solved by a computer. There are four types of unknown variables: three speed components and pressure, and the total number is 143,386 yuan.

FIG. 6 is a graph showing a convergence curve obtained as a result of a numerical experiment. The convergence determination condition is that the relative residual is 1.0e-5 or less. As the relaxation method, the SOR method (relaxation coefficient ω = 0.2) was used, and n _b = 2, n _u = 2 and n _s = 2. As a comparative example, the case where the number of iterations of the relaxation method in the AMG cycle is fixed to 4 is indicated by a broken line.

  As can be seen from FIG. 6, according to the present embodiment, convergence is started from a stage earlier than the comparative example, and after convergence has once stopped, convergence is started again after the third reset. On the other hand, in the comparative example, the convergence speed is slower than that of the present embodiment, and once the accuracy is higher than that of the present embodiment, the convergence is eventually stopped. Therefore, the present embodiment has higher convergence in a short time, and the accuracy finally obtained greatly exceeds the comparative example.

  In the present invention, the convergence can be improved in the numerical calculation using the residual cutting method adopting the AMG method as an internal solver. For example, for improving the accuracy and speed of numerical calculation such as structural analysis and fluid analysis. ,It is valid.

It is a figure which shows notionally the V-cycle as an example of an AMG cycle. It is a figure which shows an example of an AMG cycle conceptually, (a) is a W-cycle, (b) is a sawtooth-cycle. It is a flowchart which shows the numerical calculation method which concerns on one Embodiment of this invention. It is a block diagram which shows notionally the structure of the numerical calculation apparatus which concerns on one Embodiment of this invention. It is a figure which shows notionally the (nu) -cycle which is an example of the AMG cycle used in this embodiment. It is a graph which shows the convergence curve obtained as a result of the numerical experiment.

Explanation of symbols

DESCRIPTION OF SYMBOLS 1 Numerical calculation apparatus 10 1st calculating part 11 Internal solver 20 2nd calculating part

Claims (8)

When the coefficient matrix (N rows and N columns, N is a positive integer) when discretizing the partial differential equation to be satisfied by the physical quantity U is A and the inhomogeneous term (source term) is f, using a computer,
A ・ U = f
And performing a numerical calculation of the physical quantity U,
Set the initial value U ⁰ of the physical quantity U,
0 is set as the initial value for the number of repetitions m, 0 is set as the initial value of the perturbation amount φ, and (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r.
While incrementing the repeat count m,
The predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a first step of obtaining by iteration,
A modified approximate value φ ^m that minimizes the L2 norm of the residual r ^{m + 1} is obtained from the predicted approximate value Ψ ^m by the optimization routine by the second arithmetic unit, and is used as an approximate solution U ^{m + 1} (U ^m + φ ^m ) and the second step of giving (r ^m −A · φ ^m ) as the residual r ^{m + 1} is repeatedly executed until the approximate solution U ^m converges,
The first step includes
AMG method is used for the internal solver, and
In the AMG cycle, when the number of relaxation iterations n falls below the grid level, it is set to be relatively small when the level is relatively shallow and relatively large when the level is relatively deep. A numerical calculation method characterized by using a ν-cycle set to a predetermined value when going up. In claim 1,
The ν-cycle has a relaxation method iteration number n
1) For mitigation before “restricting” from level k to level k + 1:
n = max (k−n _b , 0)
2) For mitigation before “interpolating” from level k + 1 to level k:
n = n _u
(N _b and n _u are natural numbers)
A numerical calculation method characterized by that. When the coefficient matrix (N rows and N columns, N is a positive integer) when discretizing the partial differential equation to be satisfied by the physical quantity U is A and the inhomogeneous term (source term) is f, using a computer,
A ・ U = f
And performing a numerical calculation of the physical quantity U,
Set the initial value U ⁰ of the physical quantity U,
0 is set as the initial value for the number of repetitions m, 0 is set as the initial value of the perturbation amount φ, and (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r.
While incrementing the repeat count m,
The predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a first step of obtaining by iteration,
A modified approximate value φ ^m that minimizes the L2 norm of the residual r ^{m + 1} is obtained from the predicted approximate value Ψ ^m by the optimization routine by the second arithmetic unit, and is used as an approximate solution U ^{m + 1} (U ^m + φ ^m ) and the second step of giving (r ^m −A · φ ^m ) as the residual r ^{m + 1} is repeatedly executed until the approximate solution U ^m converges,
When the residual cutting rate becomes smaller than a predetermined value, the modified approximate solution φ ^m is discarded, and a reset process involving a change in the algorithm of the internal solver used in the first step is performed. Method. In claim 3,
The first step uses an AMG method as the internal solver,
In the reset process, as a change of the algorithm of the internal solver, the setting of the number of iterations n of relaxation in the AMG cycle is changed. When the coefficient matrix (N rows and N columns, N is a positive integer) when discretizing the partial differential equation to be satisfied by the physical quantity U is A and the inhomogeneous term (source term) is f, using a computer,
A ・ U = f
Is a device that performs numerical calculation of the physical quantity U,
Set the initial value U ⁰ of the physical quantity U,
0 is set as the initial value for the number of repetitions m, 0 is set as the initial value of the perturbation amount φ, and (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r.
While incrementing the repeat count m,
The predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a first step of obtaining by iteration,
A modified approximate value φ ^m that minimizes the L2 norm of the residual r ^{m + 1} is obtained from the predicted approximate value Ψ ^m by the optimization routine by the second arithmetic unit, and is used as an approximate solution U ^{m + 1} (U ^m + φ ^m ), and the second step of giving (r ^m −A · φ ^m ) as the residual r ^{m + 1} is repeatedly executed until the approximate solution U ^m converges,
The first step includes
AMG method is used for the internal solver, and
As the AMG cycle, when the number of relaxation iterations n falls below the grid level, it is set to be relatively small when the level is relatively shallow, and relatively large when the level is relatively deep, whereas the grid level A numerical computation device using a ν-cycle set to a predetermined value when going up. When the coefficient matrix (N rows and N columns, N is a positive integer) when discretizing the partial differential equation to be satisfied by the physical quantity U is A and the inhomogeneous term (source term) is f, using a computer,
A ・ U = f
Is a device that performs numerical calculation of the physical quantity U,
Set the initial value U ⁰ of the physical quantity U,
0 is set as the initial value for the number of repetitions m, 0 is set as the initial value of the perturbation amount φ, and (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r.
While incrementing the repeat count m,
The predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a first step of obtaining by iteration,
The corrected approximate value φ ^m that minimizes the L2 norm of the residual r ^{m + 1} is obtained from the predicted approximate value Ψ ^m by the optimization routine by the second arithmetic unit, and is used as an approximate solution U ^{m + 1} (U ^m + φ ^m ) and the second step of giving (r ^m −A · φ ^m ) as the residual r ^{m + 1} is repeatedly executed until the approximate solution U ^m converges,
When the residual cutting rate becomes smaller than a predetermined value, the modified approximate solution φ ^m is discarded, and a reset process involving a change in the algorithm of the internal solver used in the first step is performed. apparatus. When the coefficient matrix (N rows and N columns, N is a positive integer) when discretizing the partial differential equation to be satisfied by the physical quantity U is set to A and the inhomogeneous term (source term) is f,
A ・ U = f
And a numerical calculation program for calculating a physical quantity U,
Set the initial value U ⁰ of the physical quantity U,
0 is set as the initial value for the number of repetitions m, 0 is set as the initial value of the perturbation amount φ, and (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r.
While incrementing the repeat count m,
The predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a first step of obtaining by iteration,
A modified approximate value φ ^m that minimizes the L2 norm of the residual r ^{m + 1} is obtained from the predicted approximate value Ψ ^m by the optimization routine by the second arithmetic unit, and is used as an approximate solution U ^{m + 1} (U ^m + φ ^m ) and the second step of giving (r ^m −A · φ ^m ) as the residual r ^{m + 1} is repeatedly executed until the approximate solution U ^m converges,
The first step includes
AMG method is used for the internal solver, and
In the AMG cycle, when the number of relaxation iterations n falls below the grid level, it is set to be relatively small when the level is relatively shallow and relatively large when the level is relatively deep. A program for numerical calculation that causes a computer to execute processing using a ν-cycle set to a predetermined value. When the coefficient matrix (N rows and N columns, N is a positive integer) when discretizing the partial differential equation to be satisfied by the physical quantity U is set to A and the inhomogeneous term (source term) is f,
A ・ U = f
And a numerical calculation program for calculating a physical quantity U,
Set the initial value U ⁰ of the physical quantity U,
0 is set as the initial value for the number of repetitions m, 0 is set as the initial value of the perturbation amount φ, and (f−A · U ⁰ ) is set as the initial value r ⁰ of the residual r.
While incrementing the repeat count m,
The predicted approximate value [psi ^m of A · φ = r ^m, the first operation unit having an internal solver, a first step of obtaining by iteration,
A modified approximate value φ ^m that minimizes the L2 norm of the residual r ^{m + 1} is obtained from the predicted approximate value Ψ ^m by the optimization routine by the second arithmetic unit, and is used as an approximate solution U ^{m + 1} (U ^m + φ ^m ) and the second step of giving (r ^m −A · φ ^m ) as the residual r ^{m + 1} is repeatedly executed until the approximate solution U ^m converges,
When the residual cutting rate becomes smaller than a predetermined value, the modified approximate solution φ ^m is discarded, and a numerical value that causes the computer to execute a reset process that involves changing the algorithm of the internal solver used in the first step Calculation program.

Images