StaticLinearSolver-UG4 component for VRL by A.Vogel and M.Hoffer

StaticLinearSolverCustom.groovy

/**
 * Script to perform an Advection-Diffusion Problem
 * 
 * Author: M. Hoffer, A. Vogel
 */
import eu.mihosoft.vrl.annotation.*;

import edu.gcsc.vrl.ug.api.F_GlobalDomainRefiner;
import edu.gcsc.vrl.ug.api.I_IRefiner;
import edu.gcsc.vrl.ug.api.I_Domain;
import edu.gcsc.vrl.ug.api.I_ApproximationSpace;
import edu.gcsc.vrl.ug.api.F_Integral;
import edu.gcsc.vrl.ug.api.*;
import java.io.File;

@ComponentInfo(name="StaticLinearSolver (custom version)", category="Custom")
public class StaticLinearSolverCustom implements java.io.Serializable {
    private static final long serialVersionUID=1L;

    @MethodInfo(name="", valueName="VTU-File", valueStyle="silent", interactive = true)
    public File invoke(
        @ParamGroupInfo(group="Problem Setup")
        @ParamInfo(name="Problem-Setup", style="default")
        I_DomainDiscretization domainDisc,
        @ParamGroupInfo(group="Problem Setup")
        @ParamInfo(name="Approximation-Space", style="default")
        I_ApproximationSpace approxSpace,
        
        @ParamGroupInfo(group="Domain")
        @ParamInfo(name="Output Filename:", style="save-dialog", options="endings=[\"vtu\"]")
        java.io.File fileOut,
        @ParamGroupInfo(group="Domain")
        @ParamInfo(name="Output-Data", style="default")
        I_VTKOutput out,
        @ParamGroupInfo(group="Domain")
        @ParamInfo(name="Number Refinements", style="default", options="value=3")
        int numRefs,
        
        @ParamGroupInfo(group="Solver Setup|true|Solver Parameters")
        @ParamInfo(name="Solver", style="selection", options="value=[\"Linear Solver\", \"Conjugate Gradient\", \"BiCGStab\", \"LU\"]")
        String solverName,
        @ParamGroupInfo(group="Solver Setup")
        @ParamInfo(name="Preconditioner", style="selection", options="value=[\"Geometric MultiGrid\", \"Jacobi\", \"Gauss-Seidel\", \"Incomplete LU\", \"None\"]")
        String precondName,
        @ParamGroupInfo(group="Solver Setup")
        @ParamInfo(name="Maximal Number Iterations", style="default", options="value=100")
        int maxNumIter,
        @ParamGroupInfo(group="Solver Setup")
        @ParamInfo(name="Minimal Residuum Norm", style="default", options="value=1E-10")
        double absTol,

        @ParamGroupInfo(group="Data Evaluation")
        @ParamInfo(name="Integration Subset  ", style="ugx-subset-selection-array", options="ugx_globalTag=\"TheFile\"; minArraySize=0")
        String[] vEvalSubset
        )
        {

        logDebug("DEBUG: location 1");

        String fileNameOut = fileOut.getAbsoluteFile().getAbsolutePath();

        try {
            if (fileNameOut.toLowerCase().endsWith(".vtu")) {
                fileOut.getAbsoluteFile().delete();

                System.out.println("FILE-OUT: " + fileOut.exists());
            }
        } catch (Exception ex) {
            ex.printStackTrace();
        }

//       approxSpace.const__print_statistic();
        I_Domain dom = approxSpace.domain();
        I_IRefiner refiner = F_GlobalDomainRefiner.invoke(dom);
        for(int i = 0; i < numRefs; ++i){
          refiner.refine();
        }

        logDebug("DEBUG: location 2");
 
        if("Geometric MultiGrid".equals(precondName)){
          approxSpace.init_levels();
        }
        approxSpace.init_top_surface();
        approxSpace.const__print_statistic();

        F_OrderCuthillMcKee.invoke(approxSpace, true);

        logDebug("DEBUG: location 3");

        //--------------------------------------------------------------------------------
        //--  Algebra
        //--------------------------------------------------------------------------------
        //-- matrix and vectors
        I_MatrixOperator A = new MatrixOperator();
        I_GridFunction u = new GridFunction(approxSpace);
        I_GridFunction b = new GridFunction(approxSpace);

        logDebug("DEBUG: location 4");

        //-- create Convergence Check
        I_ConvCheck convCheck = new ConvCheck();
        convCheck.set_maximum_steps(maxNumIter);
        convCheck.set_minimum_defect(absTol);
        convCheck.set_reduction(1e-20);


        logDebug("DEBUG: location 5");

        //-- create Solver with ilu preconditioner
        I_ILinearOperatorInverse solver = null;
        if     ("BiCGStab".equals(solverName)) solver = new BiCGStab();
        else if("Conjugate Gradient".equals(solverName)) solver = new CG();
        else if("Linear Solver".equals(solverName)) solver = new LinearSolver();
        else if("LU".equals(solverName)) solver = new LU();
        else errorExit("Cannot find solver: " + solverName);
        solver.set_convergence_check(convCheck);

        logDebug("DEBUG: location 6");

        // -- create GMG Preconditioner
        if(!precondName.equals("None") && !solverName.equals("LU")) {
            I_ILinearIterator precond = null;
 
            if     ("Jacobi".equals(precondName)) precond = new Jacobi();
            else if("Gauss-Seidel".equals(precondName)) precond = new GaussSeidel();
            else if("Incomplete LU".equals(precondName)) precond = new ILU();
            else if("Geometric MultiGrid".equals(precondName)) {
                precond = new GeometricMultiGrid(approxSpace);
                ((GeometricMultiGrid) precond).set_discretization(domainDisc);
                ((GeometricMultiGrid) precond).set_base_level(0);
                ((GeometricMultiGrid) precond).set_base_solver(new LU());
                ((GeometricMultiGrid) precond).set_smoother(new ILU());
                ((GeometricMultiGrid) precond).set_cycle_type(1);
                ((GeometricMultiGrid) precond).set_num_presmooth(3);
                ((GeometricMultiGrid) precond).set_num_postsmooth(3);
            }
            else errorExit("Cannot find preconditioner: "+precondName);

            logDebug("DEBUG: location 7");
 
//            solver.set_preconditioner(precond);
            try {
                ((I_IPreconditionedLinearOperatorInverse)solver).set_preconditioner(precond);
            } catch (Exception e) {
                // in case of solver is LU
                // there will be (maybe) a CastExeption
                System.err.println(getClass().getSimpleName() +": solver is LU therefore CastException !?");
            }

            logDebug("DEBUG: location 8");
       }

       logDebug("DEBUG: location 9");

        //--------------------------------------------------------------------------------
        //--  Apply Solver
        //--------------------------------------------------------------------------------
        //-- 1. assemble matrix and rhs
        domainDisc.assemble_linear(A, b);

        //-- 2. set dirichlet values in start iterate
        u.set(0.0);
        domainDisc.adjust_solution(u);

        //-- 3. init solver for linear Operator
        F_Print.invoke("Starting Linear solver ...\n");
        if(solver.init(A) == false){
            errorExit("Initialization of Linear Solver failed.");
        }

        logDebug("DEBUG: location 10");
        
        //-- 4. apply solver
        if(solver.apply_return_defect(u,b) == false){
            errorExit("Linear Solver did not converge.");
        }
        F_Print.invoke("... Linear solver done.\n");
  
        //--------------------------------------------------------------------------------
        //--  Output
        //--------------------------------------------------------------------------------

        for(int i = 0; i < vEvalSubset.length; ++i){
            double val =  F_Integral.invoke(u, approxSpace.const__name(0), vEvalSubset[i]);
            F_Print.invoke( "Integral on Subset "+vEvalSubset[i]+":"+val+"\n");
        }

        logDebug("DEBUG: location 11");

        if (fileNameOut.endsWith(".vtu")) {
            fileNameOut = fileNameOut.substring(
                0,fileNameOut.length()-4);
        }

        logDebug("DEBUG: location 12");

        out.print(fileNameOut, u);

        logDebug("DEBUG: location 13");

        return new File(fileNameOut+".vtu");
    }

    private void errorExit(String s){
        eu.mihosoft.vrl.system.VMessage.exception(
                "Setup Error in StaticLinearSolver: ", s);
    }

    private void logDebug(String msg) {
        if(isDebug()) {
            System.err.println(msg);
        }
    }

    private boolean debug;

    public boolean isDebug() {
        return debug;
    }

    public void setDebug(boolean debug) {
        this.debug = debug;
    }
}